Genome Analysis of Three Bacteriophages Targeting Multi Drug Resistant Clinical Isolate of E. coli, Host Receptor Prediction Using Machine Learning Tools and the Evaluation of Monophages and Cocktail In Vitro and In Vivo

1. Introduction

Alternative therapy to tackle antimicrobial resistant pathogens are one of the promising ways in the current multi drug resistant era. Bacteriophages contribute to an essential part of human microbiome and around 10³⁰ phages prevail in the biosphere [1]. Owing to these, they qualify to be a tolerable, non-toxic and a reliable tool to serve as adjuvants to antibiotics or for antibiotic sparing therapy. Phage therapy did not gain therapeutic attention due to the discovery and ease of use of broad spectrum antibiotics in the 1930s – 1940s. Nevertheless, the current situation of multidrug resistance has led to resurgence and attention towards phage therapy. In 2017, WHO released a list of 12 high priority anti microbial resistant (AMR) pathogens of which carbapenem resistant ESBL producing Enterobacteriaceae was classified under the critical pathogens list. Pathogens under the critical priority list are typically curtailed by using last resort drugs. Due to evolutionary selection pressures, microbes have gained resistance to even last resort drugs like colistin through dissemination of mobile colistin resistance (mcr-1-10) plasmids both in environment and clinical settings [2,3]. Widespread dissemination of genetic elements that confers resistance to last resort drugs like colistin will catapult humanity into post antibiotic era wherein, none of the existing antimicrobials will be effective. Hence there is an urgent need to tackle multidrug resistant (MDR), XDR and Pan DR microbes using alternate approaches like phage therapy. Recently, FDA has approved the clinical trials for the use phage cocktail WRAIR-PAM-CF1 to restrain Pseudomonas aeruginosa in chronically colonized CF patients [4,5].In case of critically ill COVID patients as a compassionate measure, FDA had accorded the approval to use Phagebank ™ therapy from Advanced Phage Therapeutics to curtail MDR strains of A. baumannii, S. aureus and P.aeruginosa[6]. Phage therapy is in vogue in Georgia, Polland and Russian Federation. In Georgia, Eliavia Phage Therapy Center is successfully functioning from 2004, which provides effective solutions to patients from across the globe who do not respond to conventional antibiotic therapy[7]. In certain countries like Belgium, phages are used as Magistral therapy to treat patients when conventional therapy fails [8]. Despite these recent promising developments phage therapy is yet to gain foothold as an active therapeutic practice, due to resistance/persistence of certain strains to phage therapy, which is attributed to adaptation/species diversification [9]. A previous study reported that lytic bacteriophage belonging to Myoviridae from urban sewage, vB_Eco4M-7, was effective against multiple E. coli O157 strains, and the phage did not harbor any toxins and virulence factors [10]. Another study showed that a phage PDX, belonging to Myoviridae, isolated from wastewater in Portland killed diarrhoeagenic enteroaggregative E. coli isolates, leaving the human microbiome undisturbed [11]. Phage therapy has also been widely used as a combination therapy along with antibiotics to tackle MDR pathogens. A recent report showed that the overexpression of AmpC in E. coli promoted its susceptibility for phage lysis, possibly by employing OmpA (outer membrane protein) as receptor [12]. Attempts by microbe to enhance ampicillin resistance collaterally led to its susceptibility to phage, which favors the concurrent use of antibiotics and phage therapy in synergy. A previous study has shown that the sub lethal doses of ciprofloxacin and ECA2 phage (Podoviridae) exhibited synergy against E. coli, causing 7.8 log CFU/ml decline in 8 h [13]. In the present study, 3 different phages targeting MDR clinical isolate of E. coli U1007 was isolated from water bodies (2 Rivers and a Hospital Waste Water) and the isolated phages were characterized. From the genome sequence of phage, the Receptor Binding Protein RBP (tail spike protein) was identified and by employing machine learning algorithms, host cell surface receptor was predicted, which was further validated by growing host under conditions that is well known to alter the expression of some of these receptors. Furthermore, the ability of the isolated phage to curtail growth of XDR E. coli was explored with monophages and phage cocktail both in the presence and absence of colistin in vitro and in vivo.

2. Materials and Methods

2.1. Screening for Bacteriophages

Bacteriophages were screened against the MDR E. coli strains U3790 and U1007 strains from different sources. Genome of U1007 strain was sequenced and AMR genes were identified by Resfinder, RAST and Roary and is presented in Table S1. The antimicrobial resistance profile of these strains determined by two fold microbroth dilution assay and is provided in Table S2. The Ganges river, The Cauvery river, Cooum River, Hospital waste water, pond water samples, samples from cowsheds and soil samples from farmland were screened for phages using spot test [14]. Briefly, the samples were incubated with host culture at 37°C and after incubation for 16 – 24 h, the samples were centrifuged at 5000 rpm for 15 min. 3 - 5 µl of the supernatant harboring phages were spotted on agar plates overlaid with the host (clinical isolates of E. coli U3790/U1007) and the plates were incubated for 18 – 24 h at 37°C. Presence of bacteriophages against a specific host strain can be identified by appearance of clear zones or plaques on the agar plates.

2.2. Isolation and Purification of Bacteriophages

In order to isolate the phage, the phage containing supernatant was filtered through 0.45 µm and 0.22 µm syringe filters [15]. The filtered phage lysate was serial diluted in SM buffer (100 mM Sodium Chloride, 8 mM Magnesium sulphate, 50 mM Tris hydochloride (pH 7.5)) and each dilution was allowed to incubate with mid log cells of the host. After 20 min of incubation, the phage – host mixture was added to 5 ml of soft agar (0.7% Luria Bertani Agar) and overlaid on Nutrient Agar plates. The plates were allowed to solidify, and were incubated at 37°C for 18 – 24 h and observed for plaques. A single plaque was then picked and resuspended in 1 ml of SM buffer, serial diluted, mixed with host cells and overlaid on nutrient agar plates as mentioned earlier. The procedure was repeated three times to obtain triple purified plaques containing identical morphology. The phage titer was determined at each step and represented as PFU/ml.

2.3. TEM Imaging

The triple purified phages were enriched using the specific host (E. coli U1007 strain) to obtain a high phage titre (>10¹² PFU/ml). 10 µl of the high titer phage lysate was added to the carbon coated copper grid and was allowed to attach for 2 min [16]. The excess phage lysate was immediately removed carefully using a filter paper and then the phages were stained with 2% uranyl acetate for less than a minute. The excess stain was removed, the grid was allowed to dry and was then observed under a FEI Tecnai G² 20 S-Twin Transmission Electron Microscope (TEM) at 200 kV. The TEM images were analysed for the phage morphology to discern the family to which phage belongs.

2.4. One Step Growth Curve

Burst size and latent period was determined for the isolated phage using one-step growth curve analysis [17]. Mid log cells of the host bacteria were mixed at a multiplicity of infection (MOI) of 1 and incubated at 37°C for 5 – 10 min for adsorption. The cells with adsorbed phages were harvested by centrifugation and resuspended in Nutrient Broth and incubated at 37°C. Phage titer was determined for the samples at different time intervals until 60 min. Latent period is the time interval between the phage adsorption to the host and the host cell lysis. Burst size is the number of phages from an infected host cell and is calculated as the ratio of average PFU/ml of latent period to average PFU/ml of last three time points. The experiment was performed in triplicates and was reported as SD from the mean.

2.5. Host Specificity

The host range of the isolated phage was determined by spot assay with 10 microbes which included 8 different E.coli strains, Salmonella enterica serovar Typhimurium and Klebsiella pneumoniae [18]. 300 µl of the bacterial culture to be tested as host was added to soft agar, overlaid on nutrient agar plates and allowed to solidify. 2 µl of purified phage lysate was spotted on overlaid plates and were allowed to dry. The plates were incubated at 37°C for 18 – 24 h and then observed for clear zones. Presence of plaques represent the susceptibility of the bacterial culture to the isolated host.

2.6. Temperature and pH Stability

The temperature sensitivity of the phage was studied by incubating phage at different temperatures from 4°C, 16°C, 25°C, 37°C, 45°C, 65°C and 95°C for 1 h [19]. After incubation, the phage lysates were made to adsorb with the host and plated to determine phage titers. Similarly, pH sensitivity of the phage was analysed by incubating the phage at different pH 3.0, 5.0, 7.0, 9.0 and 11.0 for 1 h and PFU/ml was determined post incubation [20]. All experiments were performed in triplicates and represented as percent survival rate.

2.7. In Vitro Time Kill Study

In order to evaluate the efficiency of phage to inhibitthe growth of antibiotic resistant strains in vitro, time kill study [21] was performed with monophages, and different phage cocktails both with and without colistin.Mid log cells were subjected to different treatments and at regular time intervals viz., 0, 2, 4, 6 and 24 h, the samples were withdrawn, serial diluted and plated on LA plates. After overnight incubation at 37°C, the CFU/ml was calculated and the difference in colony counts were analyzed.

2.8. In Vivo Toxicity Study

The toxic effect of phage on zebrafish (Danio rerio) PET strains [22] was discerned by injecting 10 µl of purified phage lysate (10⁴, 10¹⁰, and 10¹² PFU/ml) intramuscularly and the liver and brain enzyme profiles were analysed. Phage injected fish and uninjected fish were monitored for 48 h and then the fish were euthanised and dissected. The brain and liver tissue were isolated, homogenized and the clear supernatant was used for further analysis. Brain and liver enzyme profiles were evaluated using acetylcholine iodide and α/β naphthyl acetate as substrates, respectively [23]. Any significant changes in enzyme levels relative to untreated control was deemed to be toxic to zebrafish.

2.9. In Vivo Infection

The efficacy of monophage and phage cocktail alone or in combination with colistin in preventing growth of MDR strain was evaluated using zebrafish infection [24]. 10 µl of E. coli U1007 (~10⁸ CFU/ml) was injected intramuscularly and 2 h post infection the fish were split to receive different treatment combinations viz., monophage/phage cocktail alone, monophage/phage cocktail + colistin, colistin alone. 24 h post treatment, the fish were euthanised by Immobilization by submersion in ice cold water, the infected muscle tissue was dissected, homogenized, serially diluted, and plated on LA plates. After incubation for 24-48h at 37^oC, colony counts were determined and represented as mean CFU from triplicate values.

2.10. Genome Sequence and Analysis

The phage DNA was extracted using Cetyl Trimethyl Ammonium Bromide (CTAB) DNA precipitation method [25]. 1.5 ml of purified high titer phage lysate was incubated at 22^oC for 15 min with 10 ng of RNase A and 10 U of DNase I. Post incubation, 80 µl of 0.5 M EDTA and 50 µg of Proteinase K was added and was maintained at 45^oC for 15 min, followed by addition of 5% CTAB and incubation in ice for 15 min. The precipitated DNA was harvested at 8000 g and the pellet was resuspended in 1.2 M NaCl. The DNA was further precipitated and washed with ethanol. The pellet was air dried and suspended in 10 mM Tris buffer, which was stored at -20^oC until use. The DNA was quantified using Qubit Fluorometer and was then sequenced using Oxford MinIon Nanopore sequencer. Library preparation was performed as per manufacturer’s instruction using Ligation sequencing kit (SQK-LSK109) and the sequencing run was performed without live basecalling. Basecalling and demultiplexing were performed using Guppy and de novo assembly of the reads were done using Canu Assembly software for the U1G genome [26]. For the CR and M phage genomes, basecalling and demultiplexing of the Minion raw data were done using Guppy. Adapters were removed using Porechop (https://github.com/rrwick/Porechop) and further processing using fastp[27] and only reads above 200 bp were retained. De novo assembly of the reads was done using Flye (https://github.com/fenderglass/Flye) with asm_coverage of 50 and four polishing iterations using an expected genome size of 60 Kb. The genome assembly and completeness were assessed using QUAST Version 5.2.0(https://github.com/ablab/quast). The assembled reads were annotated using Rapid Annotation using Subsystem Technology (RAST) [28]. The assembled genome was searched against nucleotide database available at National Center for Biotechnology Information (NCBI) using Basic local alignment search tool (BLASTn) to identify close homologs and the top phage hits with query coverage above 80 were considered for the analysis [29]. Presence or absence of tRNAs were identified using tRNAscan-SE search server [30] and the reads were also fed to PHASTER to determine intactness [31]. The Average Nucleotide Identity (ANI) between the close homologs of U1G was calculated using the ANI/AAI-Matrix tool (http://enve-omics.ce.gatech.edu/g-matrix/) [32]. Phylogenetic tree of phoage homologs was obtained using ANI-Distance clustering method, UPGMA. For RBP homologs, tree was constructed using the Maximum Likelihood method. The tree was visualised using iTOL v6(https://itol.embl.de/) [33]. Antibiotic resistant and virulent genes in the phage genome were analysed using ResFinder and Virulence Finder of CGE [34,35]. Potential RBP sequences in the annotated genomes of all three phages were identified by conducting a BlastP search of the amino acid sequences of the assembled bacteriophage genomes against the RBP protein sequence database locally [36].

2.11. E. coli Bacteriophage Host Receptor Prediction Using Machine Learning Based Multi Class Classifier

The initial entries of the RBP nucleotide and protein sequences were obtained from the database curated by a study [36] conducted to predict bacterial hosts which consisted of 1232 RBP sequences belonging to nine different bacterial hosts, out of which 400 had E.coli as the bacterial host. The bacteriophages that targets E.coli as the host were mapped to the PhRED database [37] which provides information on the corresponding receptor proteins involved in the bacteriophage host interaction of which only 71 RBP entries had the receptor information. For RBP database entries with PhRED receptor information, the sequence information was updated by a manual search of annotated proteins and CDS from the genome sequence of bacteriophages. Common search terms used were ‘tail fiber protein’, ’tail protein’, ‘tail fiber’ along with the bacteriophage name. RBP database entries with sequence information pertaining to the receptorinformation were updated from a literature survey of experimental studies conducted on bacterial host receptor interactions [38,39]. The curated database of RBP and receptor information with 160 entries were used for building the ML model (Table S3). The RBP sequences were presented as a vector of numerical features extracted from both the nucleotide and protein sequences based on the script available with the study [36]. A total vector of 218 numerical features (Table S3) was retrieved for each of the RBP dataset entries.

Data preprocessing involves checking for null values in the dataset, obtaining a balanced dataset, and training it. Upon performing exploratory data analysis, it was found that the RBP dataset is imbalanced in nature in terms of the target variable distribution with some labels having just a single entry. In order to address this issue, only RBP entries with label count two or more were included and RBP entries with single-entry outliers were removed, setting the final size of curated dataset to 155. We included SMOTE (Synthetic Minority Oversampling Technique) [40], a method of oversampling that creates artificial samples from the class with the lowest count. For training the classifier, SMOTE is utilized to create a training set that artificially balances the class distribution. Three classification algorithms, Random Forest, Multinomial Logistic Regression and Decision Tree were used for the construction of the multiclass classification models and were validated using the Nested Cross Validation approach. The range of parameters to be run for the classification models was defined. Multinomial Logistic Regression: C: [0.001, 0.01, 0.1,1,10,100,1000]andDecisiontreeClassifier: max_features”:[2,4,6],”criterion”:[“gini”,”entropy”],RandomForestClassifier: n_estimators”:[10,50,100], max_features: [5,10,15], criterion:[“gini”,”entropy”] in order to run the Nested Cross Validation with an Outer Loop Fold of 10 and Inner Loop Fold of 5. In this study, we implemented two feature selection methods, Analysis of Variance (ANOVA) and L2 Regularization, and compared the performance of the trained model with those two datasets to that of the 218 Features incorporated dataset.

The performance evaluation scores like accuracy, precision, F-1 Score and Matthew’s Correlation coefficient (MCC) were calculated from the three categories of datasets namely with all features selected dataset with 218 features and ANOVA feature selected dataset with 30 highest scoring features, and L2 -Regularization selected Dataset with 110 features trained with the three classification algorithms. The ROC Curves were generated for each Dataset - Classification algorithm combination. The implementation of Classification algorithms and cross validation was done using the Scikit-learn package (version 1.0.2) available in Python[41]. After optimizing the ML model, the RBP sequences of the phages were used to predict the potential host receptors of the bacteriophages.

3. Results

3.1. Identification of Phage against XDR E. coli Strains

Different sources viz., The Ganges River, The Cauvery River, Cooum River, Hospital waste water, pond samples, soil samples from farmland, samples from cow shed were collected and screened for presence of bacteriophages specific to the multidrug resistant clinical isolate of E. coli ( U1007). Genome sequencing of U1007 E. coli strain (SRA Accession # PRJNA988283) followed by the identification of antimicrobial resistant genes (ARG) revealed that the genome of U1007 harbored resistance to 48 antimicrobials belonging to multiple antimicrobial classes viz., cephalosporins, floroquinolones, macrolides, carbapenems, tetracyclines including efflux transporters that have been associated with colistin resistance (Table S1). Determination of MIC using two fold microbroth dilution assay also showed that strain showed resistance to multiple antimicrobials belonging to 7-8 different antimicrobial classes and displayed an MIC of 4µg/ml for colistin (Table S2). The drug resistance profile of these strains were also reported earlier [42]. Diverse water samples incubated with the host bacterium (U1007 and U3790) were centrifuged and the supernatant was spotted on plate overlaid with 10 different bacteria (8E. coli strains and Salmonella enterica and Klebsiella pneumoniae). We found that the water samples from The Ganges River (U1G) The Cooum River (CR) and Hospital waste water (M Phage) harbored lytic bacteriophages specific to the XDR strain U1007, U1G showed faint lysis against U3790 strain whereas the other two phages (CR & M) displayed lytic potential against another clinical isolate U2354 (Figure S1 & Figure S6).

3.2. Isolation and Purification of U1G, CR and M Phages

Spot assay revealed the presence of bacteriophages targeting U1007 strain from The Ganges river e, Cuoom River and Hospital waste water. In order to isolate these phages, the supernatant was filtered, serial diluted, incubated with U1007 for 20 min and overlaid on nutrient agar plates. Post incubation, the plaques were identified based on clear zone of lysis. The phage morphology was carefully observed and a single plaque morphology was taken up for purification. The phages were triple purified ensuring that similar plaque morphology was repeatedly obtained (Figure 1) and a phage stock of high titer (10¹³ PFU/mL) was stored at 4⁰C for further use. The phage specific to U1007 from Ganges was designated as U1G, Cuoom river as CR and Hospital waste water as M phage respectively.

3.3. TEM Imaging

The morphology of phages were observed by staining the phages with 2% uranyl acetate using Transmission electron microscopy. The TEM image showed that all 3 phages possessed an icosahedral head of varying diameter ranging from 71.76 nm for U1G, 31.64 nmfor CR phage and 54.03 nm for the M phage. U1Ghad a very short non-contractile tail, whereas CR phage had a non- contractile tail of 24.37 nm and M phage had a long non-contractile tail of 104.13 nm (Figure 1). Thus based on TEM morphology U1G belongs to Podoviridae, CR phage belongs to Myoviridae and M phage belongs to Siphoviridae.

3.4. Genome Analysis

Whole genome sequencing of U1G, CR and M phages were performed using Oxford MinIon Nanopore platform. The raw reads of the U1G phage were trimmed and corrected using Porechop and assembled to contigs using Canu v 1.8. The assembled contigs were annotated using RAST and the genome sequence was submitted in NCBI GenBank (Accession Number: MZ394712). The genome of U1G is 73,275 bp long and it has a GC content of 43%, N50 of 73275 and L50 as 1. RAST annotation revealed the presence of 91 coding sequences, out of which 31 sequences code for phage packaging, replication, and other functions (Table S3). Rest of the proteins were annotated as hypothetical. Annotation of phage genome sequence by PHASTER revealed that the genome was intact. PhagePromoter tool predicted 54 promoters in the genome and absence of tRNAs were revealed by tRNAscan-SE. Genome map of U1G constructed using SnapGene 6.2 showed that the majority of the identified functional modules clustered in the first half of the genome (Figure 2A).

The assembled contigs of CRphage and Mphage were annotated using RAST and the genome sequence was submitted to NCBI GenBank (Submission Ids 2708670 & 2708284 Accession Numbers: OR061068 & OR061069). The assembled genome of CRphage formed a circular genome with 45,236 bp with coverage of 248X. The maximum coverage for M phage was 21X. From the Flye assembly of M phage reads, the contig with maximum coverage alone considered for further analysis. The selected contig was a circular one with 45,294 bp long. The GC content of CR and M phages was found to be 50.9%. RAST annotation identified the presence of 95 and 104 features in CR and M phages, respectively(Tables S4 and S5) and both the genomes were marked by the absence of tRNAs as revealed by tRNAscan-SE. Genome map of CR and M phage is reported in Figures 2B and 2C, respectively.

The top 14 close homologs of U1G phage genome from BLASTn search were, PGN829, vB_EcoS_Uz-1, vB_EcoM_PD205, PD38, Bp4, St11Ph5, vB_EcoP_PhAPEC7, vB_EcoP_PhAPEC5, vB_EcoP_PhAPEC7, vB_EcoP-ZQ2, Caudoviricetes sp. Isolate 355, vB_EcoP_G7C, vB_Eco_F22, ECBP1, and E20. The CR and M phages showed 97% similarity and there are seven close homologs identified for both phages and are Caudoviricetes sp. Isolate ctNve1, Caudoviricetes sp. Isolate ctMtm1, 0116_121510, Caudoviricetes sp. isolate ctkLI7, vB_EcoS_SCS92, NTEC3, and Bacteriophage sp. isolate 2894_61690. The close homologs of the three phages were used for the ANI matrix calculation and the construction of Phylogram (Figure 3). The results indicate that the U1G phage is distinct from CR and M phages as they form two separate clusters. NCBI BLASTn results showed that the U1G genome matched 95.98% with Escherichia phage PGN829.1 with 90% query cover and e value 0.

Phylogenetic tree also revealed that U1G is a close relative to Escherichia phage PGN829.1 (Figure 3B). However, PGN829.1 was classified under a new family Schitoviridae [43], which is highly similar to Podoviridae in morphology but harbors virion associated RNA polymerase. Mauve based comparative genome analysis of the close homologs of U1G (identified using ANI matrix) revealed three homologous blocks shared among the genomes. The homologous blocks are syntenic but positional variation for these syntenic regions was observed in different phages (Figure S2A) implying that these phages have a common ancestry. Genome alignment of CR and M phages along with their three close homologs revealed the presence of four homologs blocks shared among the genomes(Figure S2B). PHASTER showed the presence of RNA polymerase subunit in the

U1G genome and NCBI blast to detect the presence of virion associated RNA polymerase gene with U1G revealed the presence of homologous region (45314 bp – 48093 bp). Hence U1G was also classified as Schitoviridae although the presence of RNA polymerase in the U1G genome was not detected by RAST annotation. The CR and M phages showed highest similarity with Caudoviricetes sp. isolate ctkLI7 (Figure 3).

The phage tail spike protein with 107 amino acid sequence is recognised as the RBP encoded by U1G phage which is used by the phage to attach to bacterial cell surface and is highlighted (with red colour) in Figure 2A. The tail spike protein of U1G phage showed 85% similarity with the tail fibers protein of Salmonella phage SP1 of RBP database. The homologs of the tailspike protein of U1G from other bacteriophages were obtained using web BLASTp search and MEGA11 was used to construct the phylogenetic tree of tailspike protein using UPGMA method with 1000 bootstrap values. The results revealed that tailspike protein of U1G was a close homolog of PGN829 hypothetical protein and is distantly related to tail fiber protein of Salmonella phage SPHG3 (Figure S3A). The hypothetical protein of PGN829 with Accession no. AXY82585 is a lengthy one with 628 amino acids. The tail fiber domain-containing protein of Shigella phage pSb-1 (119 amino acids) and the tail spike protein of Dompiswa phage TSP7_1 (102 amino acids) are the RBPs with length comparable to that of tailspike protein of U1G.

The Potential RBP sequences of CR and M phages were identified by using local BlastP search of the RAST annotated protein sequences in the phage genome with the RBP database [36]. There were two potential RBPs identified from the CR phage and are highlighted in the genome map in Figure 2B. The first one is having 828 amino acids and is annotated as phage protein which showed 99% similarity with putative tail protein of Escherichia phage ST2 in the RBP database. The second one is with 481 amino acids and is annotated as phage tailspike protein by RAST which showed 81% similarity with the tailspike protein of Salmonella phage FSL SP-049. The tail protein from the Caudoviricetes sp. identified as the close homolog of the first RBP from the CR phage using web BlastP search(Figure S3B). Tailspike protein from Caudoviricetes sp. and Escherichia coli phage are identified as the close homologous sequences of second RBP sequence of CR phage from Blastp search (Figure S3B).

There were three potential RBPs identified from the M phage which includes a phage protein, phage tailspike protein and a hypothetical protein (Figure 2C) with lengths 852, 110 and 437 aa respectively by search against the RBP database. Putative tail protein of Escherichia phage ST2, tailspike protein of Salmonella phage FSL SP-049 and tail fiber protein_Acinetobacter phage Petty from RBP database showed 99%, 86% and 45% similarities with the three potential RBPs from M phage. Tail protein from Escherichia phage ST20 identified as the close homolog of the first RBP of M phage using web BlastP search. The tailspike protein from Caudoviricetes sp. identified as the close homolog of the second and third RBPs from M phage which are the tail spike protein and the hypothetical protein (Figure S3C).

3.5. A multiclass-classification model for bacteriophage host-receptor prediction in E. coli:

ML based multi-class classifier using phage RBP sequence was created to predict the potential host receptors for Escherichia phage U1G. The RBP host receptor dataset for E. coli targeting phage with 160 entries were obtained from the following sources viz., report by Boeckaerts, D. et al.[36], PhReD database and literature survey (Table S2). The 218 nucleotide and protein sequence features of the collected RBP sequences were obtained using the script reported earlier [36] and were used for training the ML classification algorithms. The host receptors available in the dataset includes, LPS, Tsx, OmpC, OmpA, LPS core, FhuA, LPS O antigen, LamB, OmpF, FadL, TonB, BtuB, and pili tips.The distribution details of these receptors are reported in Figure S4. Three Classification algorithms, Random Forest, Multinomial Logistic Regression and Decision Tree were used for the construction of the multiclass classification models which can predict the host receptors of E. coli targeted by phage based on the Receptor Binding Protein (RBP) sequence of the bacteriophage. Performance comparison of the classification algorithms is reported in Table 1 and the model using the Random Forest Classifier Algorithm reported the best Performance in terms of all the metrics followed by the Multinomial Logistic regression Model. We have applied two feature selection methods, ANOVA and L2 Regularization and the top selected 30 features from the first and 110 features from the second method wereused for the construction of feature selected ML models and the results are included in Table 1. Overall, the ML model using Random Forest Classifier on the data set consisting of all the 218 features selected produced the highest performance of 93% in terms of precision, 90% accuracy and an average AUC value of 0.99 in terms of individual class contributions aggregate.

The predicted RBPs from the U1G, CR amd M phages showing high sequence similarity with the RBP database entries were used for predicting the host receptor using the multiclass-classification model constructed. The tailspike protein from U1G phage, phage protein from CR phage and phage protein from M phage were used for the host receptor prediction. The host receptor prediction results from the Random Forest multi-class classifier using all features and selected features from ANOVA and L2 Regularization methods are reported in Table 2. For the U1G phage, LPS O antigen, OmpC or LamB were predicted as the host receptors. FhuA or OmpC were the predicted receptors for CR phage and OmpC, LamB, FhuA, TonB or OmpF were the potential host receptors for M phage (Table 2).

The tailspike protein was identified as the RBP protein for the Escherichia phage U1G, CR Phage and M Phage and the amino acid sequence of the tail spike protein was used for predicting the host receptor using the multiclass-classification model constructed. For CR and M phages additional RBPs identified, Phage protein in case of CR phage and Phage and hypothetical protein for M phage. All the predicted RBPs were used for host receptor prediction. The host receptor prediction results from the multi-class classifier are reported in Table 2. Two of the three ML tools predicted LPS O antigen as the host cell surface receptor when all features were selected for analysis. L2 regularization resulted in identification of OmpC/LamB/pili tips. The DecisionTree algorithm identified FhuA and TonB as potential host cell surface receptors. In case of the CR phage, the phage protein with 828 aa that had the highest similarity to the RBP and ML tool identified FhuA receptor with RF model including all features. For the M phage, the phage protein with 852 aa showing highest similarity to the RBP and identified OmpC/LamB as the probable host cell surface receptors using RF model with all features. (Table 2).

• Physiological Validation of host OmpC as the receptor for U1G phage

As RF algorithm using L2 regularization identified OmpC as one of the plausible host cell entry receptor for U1G, we checked whether differential expression of OmpC will affect U1G plaque titers, when the phage to host bacterial ratio was maintained constant. Based on the earlier well established reports,[44,45] on Osmolality inducing differential expression of OmpC, we increased the medium osmolaltity and evaluated impact of increased OmpC overexpression on the phage titers relative to medium with low osmolality. As expected, U1G phages at 10⁴ dilution when mixed with 0.4 OD (A600) of U1007 E. coli strain grown in low osmolality medium had countable plaques whereas same dilution of phages (10⁴) when mixed with 0.4 OD (A600) of U1007 E. coli grown in high osmolality medium completely lysed the bacteria resulting in complete clearance proving high phage titers (Fig 4). Thus, increased OmpC expression (induced by high osmolality) led to increased phage adsorption and cell lysis relative to cells exhibiting lower OmpC expression (caused by growth under low osmolality) which implies that OmpC as predicted by the ML tools is a potential host cell surface receptor employed by U1G.This observation would be further validated using OmpC knock out strain in an isogenic background in future studies.

3.6. Burst Size, Latent Period and Host Specificity

One-step growth curve showed that U1G possessed a relatively short latent period of 20 min, a rise period of 20 min and a burst size of 124 PFU/cell. M phage also had a short latent period of 20 min and a burst size of 150 PFU/cell.) Whereas CR phage had a short latent phase of 15 min and a moderate burst size of 117 PFU/cell (Figure S5). The host range of U1G/CR/M phages were determined by testing the lytic activity of phage against different E. coli clinical isolates U3790, U1007, U3176, IDH09733, U2354, U1024, IDH09519, and MG1655, Klebsiella pneumoniae and Salmonella Typhimurium (since the tail spike protein of Salmonella phage showed higher sequence similarity with the tail spike protein of U1G phage)by spot assay. The results showed that all 3 phages U1G/CR/M showed tropism primarily towards U1007 and mild lytic activity against U2354 but not against other isolates (Figure S6). A faint lysis zone was exhibited by U1G against U3790, though it was not significant. CR and M phages also exhibited varying ability to show lytic potential against another clinical isolate of E. coli (U2354).

3.7. pH and Temperature Sensitivity

Temperature and pH stability of U1G /CR/M phages were studied by incubating at different temperatures respectively. U1G remained stable in the temperature range of 4°C to 45°C and at 65°C, 50% loss in viability was observed (Figure 4). Whereas the other two phages CR and M were stable only until 45°C and both of them lost their viability when incubated at 65°C with CR phage being more sensitive to 65°C than M phage (Figure S7). The results of pH sensitivity revealed that U1G was stable at pH 5.0, 7.0 and 9.0 (Figure S8)at extreme pH of 3.0 and 11.0, U1G lost its viability. Altered ompC expression at extremes of pH might account for negligible lytic activity exhibited by U1G phage. CR phage was also inactive at pH 3.0 but relative to U1G phage, it retained 25 % stability at pH 11.0 and M phage was stable until pH 9.0 and at pH 11.0 it retained 25 % viability. Thus U1G had better thermotolerance whereas CR and M phages had better tolerance at alkaline pH. Thus all 3 phages differed in terms of their pH and Temperature stability which reinstates that although they target the same host, the three phages vary in their characteristics.

3.8. Time Kill Study

Time-kill was performed to evaluate the efficacy of Phage cocktails against U1007 individually and in combination with antibiotics. Our observations with monophages (U1G/CR/M) revealed that in a time kill experiment, monophages caused a decline in cell count only for the initial 2-3 h beyond which almost a complete regrowth was observed by 24h (Figure 5A) which could be attributed to the probable clonal expansion of the ressitant mutant. Treatment with phage combinations displayed a time kill trend similar to monophages but importantly, cell counts plateaued around 6h and did not increase further even after 24h which shows that phage cocktails are highly effective in restricting regrowth relative to monophages. By 24h, a 3 log decline in cell counts were observed for phage combinations(U+M+ Col) relative to the untreated control (Figure 5B) which signify that phage cocktails are indeed effective in restricting bacterial growth. Phage combinations with colistin had a slightly higher decline in cell counts than phage cocktails without colistin, which implies that despite being effective, the tested phage cocktail might have additive effect but not synergistic effect with colistin in restricting U1007. )

As the strain displayed decline in cell counts due to colistin treatment in the time kill assay (Figure 5), which was resembling trend reported earlier for colistin heteroresistant bacteria[46], we evaluated whether the MDR clinical isolate of E. coli U1007 displays colistin heteroresistance by performing population analysis profile as reported earlier [47]. U1007 strain was serially diluted and spotted on plates containing increasing concentrations of colistin and cell counts were compared with strain plated on antibiotic free plates and if the ratio greater than 0.0001 the tested strain is deemed as colistin heteroresistant. Our calculations showed (Figure S11) that the ratio of cells grown on colistin (4µg/ml) containing plates relative to colistin free plates is 0.000492 hence MDR U1007 strain in the present study was indeed colistin heteroresistant.

3.9. In Vivo Toxicity and Infection

As a representative of 3 phages, U1G phage wasevaluated for its toxicity in zebrafish by estimating the brain and liver enzyme profiles. Different titers of phages (10⁴, 10¹⁰ and 10¹² PFU/ml) were injected intramuscularly and the enzyme profiles were estimated. It was observed that there was a slight increase in α-naphthol release corresponding to control, which was not statistically significant (Figure S9A). The β naphthol (Figure S9B) release and amount of acetylcholine esterase (Figure S9C) was similar to that of untreated control. Hence the phages are unlikely to pose any toxicity to zebrafish even at a relatively high doses.

In vivo infection study was performed with monophages/phage combinations either alone or with sub MIC levels of colistin to restrict U1007 infection. The results showed that treatment with monophages resulted only upto 1.3 to 1.5 log reduction in bioburden whereas monophage +colistin treatments caused a maximum of upto ~2.1 log CFU decline relative to the untreated control. Colistin treatment alone caused 1.2 log decline in CFU (Figure S9). Phage combinations along with colistin especially U1G+ M+ Col resulted in 3.3 log decline in CFU and U1G+CR+M+Col resulted in 3 log decline in cell counts (Figure 6) which is statistically significant. Thus both bi-phage (U+M) and Triphage combination (U1G+M+CR) along with colistin caused significant 3 log decline in cell counts (Figure 6) which reiterates potential of phage cocktail and antibiotic to restrict growth of clinical isolate of E. coli in vivo. In the bi-phage (U+M) combination colsitin exhibited synergy with the phages as bi-phage (U1G+M)alone showed only 0.8 log reduction in CFU, whereas bi-phage + colistin caused ~3.3 log reduction in CFU (Figure 6). On the otherhand colistin + triple phage (U1G+CR+M) combination resulted only in a modest 0.5 log reduction in CFU. Hence ability of antibiotic to synergise with the phage is dependent on phage combination being employed.The results show that Phage cocktail along with colsitin is quite effective in curtailing bioburden of colistin heteroresistant MDR clinical isolate of E.coli in zebrafish infection model.

4. Discussion

The “antibiotic pipeline” had a tremendous growth in the late 20^th century and since the beginning of 21^st century, the pipeline is experiencing a long lag phase owing to the evolution of antimicrobial resistance leaving drugs like tigecycline, carbapenems and colistin as last resort drugs to treat infections caused by MDR pathogens [48]. Resistance to these drugs have also evolved and hence search for an alternative or adjunct therapy is the need of hour. Desperate need to tackle this situation has reminded the researchers of a century-old bacteriophage therapy. Felix d’ Herelle, who identified and termed bacteriophages, was the first to use phage to clinically treat bacterial dysentery in 1919 [49]. Another trial phage therapy was carried out by him on the people of Punjab, India to treat cholera and found a 90% reduction in mortality [50]. Nevertheless, owing to various controversies and rapid progress in the development of antibiotics pushed phage therapy behind. However, rising antibiotic resistance has re-kindled the interests towards bacteriophages. Bacteriophages have been reported against different pathogenic strains of Pseudomonas aeruginosa, Clostridium difficile, Vibrio parahemolyticus, Staphylococcus aureus, Acinetobacter baumannii, E. coli and Klebsiella pneumoniae, individually or in combination with antibiotics, thus favoring reuse of antibiotics [48]. In this study, we have isolated phages targeting MDR E. coli clinical isolate (U1007) that is resistant to multiple antibiotics belonging to 7 different classes (Table S2), including carbapenems and is ESBL positive (Figure S11B) and hence falls under the critical priority pathogen list as designated by WHO [51]. Attempts to identify phages against another colistin resistant E. coli U3790 was unsuccessful owing to the presence of the capsule as reported earlier and also due to the intact prophages within its genome [42]. Prophages in bacterial genome can evolve mechanisms like blocking phage genome injection, blocking phage binding and preventing the interaction of phage receptor on the bacterial membrane to evade superinfection by other related phages, though the exact mechanism is still not known [52]. Nevertheless, bacteriophages targeting MDR E. coli U1007 were isolated from The Ganges River (designated as U1G), Cooum River (CR) and Hospital Waste water(M) (Figure S1). We recently reported a phage KpG, belonging to Podoviridae, specific to MDR K. pneumoniae from Ganges, which was able to curtail the host’s planktonic and biofilm mode of growth [53]. There are other numerous reports available on rich diversity of bacteriophages against various pathogens being isolated from The Ganges [54,55]. This is attributed to the origin of Ganges The Himalayan permafrost, which has trapped bacteriophages from a long period and is released gradually while melting and hence it forms a seed source of bacteriophage [56]. Both Cooum River and Hospital wastewater are likely to harbor a lot of MDR bacteria. Depending upon the origin hospital waste water is likely to harbor MDR microbes from 0.58 to 40 % [57] As the prevalence of drug resistant microbes are likely to be higher in hospital wastewater, the propensity to harbor phages that target drug resistant microbes are also high. TEM imaging revealed that U1G belongs to the Podoviridae family, containing an icosahedral head of mean diameter 71.68 nm and a short non-contractile tail. Whereas CR phage possessed a shorter head of 31.64 nm and a moderate tail of 24.37 nm and is likely to belong to Myoviridae. On the other hand M phage had a head of 54 nm and a long tail of 104 nm (Figure 1) hence it probably belongs to Siphoviridae. Thus in this study we had successfully isolated phages belonging to three different families yet targeting the same host. A recent study has shown that nettle manure harbored phages belonging to Siphoviridae and Podoviridae targeting the same plant pathogen Pseudomonas syringae pv. tomato [58]. However, genome analysis of U1G by PHASTER revealed the presence of RNA polymerase and hence U1G is deemed under Schitoviridae. Reports show that by 2020 there were only 115 members that were classified under this newly proposed family Schitoviridae [43]. Interestingly we also observed during whole genome BLAST analysis (Table S4) that a prophage in Enterococcus faecium strain ME3 chromosome displayed 95.94 % identity (with 91% coverage and e-value of 0) with U1G phage genome but the presence and probable expression of antirepressor protein in U1G (Figure 2) might favor its lytic life cycle. Presence of a highly homologous prophage genome in Enterococci imply that U1G might possibly use Enterococcus faecium as a host, we tried to infect a reference strain of Enterococcus faecium with U1G using spot test, our attempts to infect was unsuccessful probably because genome harbored prophage prevents super infection by a similar phage. One step growth data showed that U1G possessed a burst size of 124 PFU/cell and a latent period of 25 min. M phage also exhibited a shorter latent period of 20 min but its burst size was 150 PFU/cell Whereas CR phage had a short burst size of 114 PFU/cell (Figure S5). The exponential growth of phage and its lytic efficiency is majorly dependent on latent period and burst size. Larger burst size usually results from a long latent period and vice versa. However, a phage with shorter latent period and relatively good burst size possesses an enhanced capability to lyse host cells faster [59]. Thus among the 3 phages employed M phage has both a short latent period and a relatively high burst size (Figure S5). Experimental evidences and predictive modeling indicate that host cell densities and phage latent periods are inversely related and phages could be evolving towards a latent period optimum which tends to maximize the population of phages that grows in the presence of a specific quality and quantity of host cells [59]. Host specificity assay revealed that all three phages U1G, CR and M were highly specific to its host U1007 and U1G displayed a faint lysis zone against U3790 whereas CR and M phage displayed faint lysis against U2354 strain (Figure S6). No clearance or appearance of plaques were seen for clinical isolates, including the reference strain MG1655. As U1G exhibited a faint zone against another colistin resistant strain U3790, this led to a hypothesis that the phage U1G might be specific to colistin resistant E. coli and does not affect colistin sensitive strains of E. coli which includes both pathogens or commensals. As colistin resistance confers chemical alteration to LPS [60], it is likely that altered LPS could serve as a receptor for phage entry, which was also predicted by ML algorithm (Table 2). Previous report showed that a panel of colistin resistant K. pneumoniae were more susceptible to lytic phage, isolated from sewage water, than to their respective colistin susceptible strains [61]. The isolated phages were observed to be negatively charged and since colistin resistant strains possess reduced cell surface negative charge, electrostatic interaction might have favored the enhanced susceptibility of colistin resistant strains to the phages[61]. In the current study, we employed a machine learning-based approaches for host receptor prediction, specifically based on the RBP sequence. To the best of our knowledge, this is the first reported study on host receptor prediction utilizing this ML-based approach. Most studies primarily focus on predicting the hosts of bacteriophages rather than the receptors [36,62]; The advantage of using machine learning tools to predict host cell surface receptors is that it will reduce significant time and labor as conventional phage adsorption studies to identify host receptor will involve laborious screens with large mutant libraries and for clinical isolates such mutant libraries needs to be created first which is cumbersome [63]. Machine learning tools will considerably reduce the labor by pinning down on a handful of putative receptors for which knock outs (even in clinical isolates) can potentially be developed and screened. Random forest algorithm using all sequence features and features based on ANOVA approach identified LPS O antigen as the host receptor employed by the phage tail protein for entry into the host. As LPS O antigen knock outs are difficult to generate in an isogenic background, we attempted to physiologically validate other receptor(s) that are predicted by the other feature selection approach. Random forest algorithm (upon L2 regularization) identified OmpC as one of the host cell surface receptors for U1G phage. As OmpC expression is modulated by medium osmolality [44,45], to validate the role of OmpC as a host cell surface receptor, U1G phage was exposed to cells grown under conditions that either upregulate OmpC (high osmolality) or downregulate OmpC (low osmolality) and equal ratio of phages to bacterial cells were maintained for both treatments and the results showed that phage titers were indeed high when OmpC is upregulated (Figure 4) this validates the prediction by RF algorithm (Table 2) that OmpC is the host cell surface receptor for the U1G phage. Future studies will attempt to create knock outs of predicted receptors (LPS O antigen, LamB, FhuA, OmpC, OmpF, TonB) in isogenic (U1007) background and reaffirm predictions of host cell surface receptors made by ML algorithms for all the three phages. The current ML algorithm takes into account only the nucleotide and protein sequence features of RBPs. However, further enhancements can be made by incorporating structural features once the relevant data becomes available[64]. Relative to CR and M phages, U1G phage was more thermostable and it retained 50 % viability at 65°C (Fig S7). Conversely M and CR phage displayed greater viability in alkaline pH relative to U1G phage (Fig S8). Interestingly, pH is also known to affect the expression of OmpC [65]. A previous report has shown that acidic pH stimulates OmpC whereas alkaline pH stimulates OmpF [64] as OmpC is the receptor employed by U1G, reduced phage titers at alkaline pH for U1G might as well be attributed to reduced expression of OmpC at alkaline pH.

In vitro time kill study with Phage combinations, individually and in combination with colistin showed that phage combination along with colistin showed better ability in restricting regrowth relative to treatment with monophages (Figure 5). Usually one would expect a drastic reduction when using phage combinations, as all three phages belonged to different families as evident from the morphology (Fig 1). But phylogenetic analysis based on genome similarity showed that U1G is distinct from CR and M which are quite closely related to each other (Figure 3). Earlier studies have shown that co-infection by different phages on the same host might result in smaller burst size and infection exclusion [66] which might account for the relatively lower titers observed in triple phage combinations.LPS O antigen is predicted as one of the targets for the phages by ML tools (Table 2), alteration of O antigen is a common phenomenon so bacteria gaining resistance to U1G is possible which necessitates the use of phage cocktails. It is likely that the bacteria during its attempt to develop phage resistance might partially lose resistance to colistin and hence the combination can achieve enhanced killing than when individually treated. Zebrafish model has been used to study the efficiency of bacteriophages in curtailing infections caused by P. aeruginosa, K. pneumoniae, E. coli and E. faecalis[53,67,68,69]. Nevertheless, the majority of studies have compared the effect of antibiotics and phage therapy and very few reports have studied the combination of antibiotic and phage therapy. In our earlier study, we found that the combination of Streptomycin and KpG (Podoviridae phage specific to K. pneumoniae) curtailed the infection by 98% relative to untreated control, whereas KpG alone caused 77% reduction and only streptomycin resulted in 63% reduction in colony counts [53]. In the present study, a drastic improvement in phage titers of upto 2.2 log CFU was observed when bi-phage (U+M) combination was used along with colistin relative to bi-phage treatment alone in fish infection study (Fig 6). A similar trend was not observed when all three phages were used in combination with colistin, which resulted in modest 0.5 log CFU difference between phage cocktails with and without colistin (Figure 6).Enhanced phage activity in the presence of subinhibitory concentrations of antibiotic termed as phage antibiotic synergy (PAS) as reported by Comeau et al., 2007 [70] was observed by many others as reviewed by North et al.,2019 [71] and is attributed to enhanced burst size or collateral sensitivity to antibiotics due to phage resistance [72]. Future studies of one step growth curve with sub MIC levels of colistin can unravel whether enhanced burst size triggered by the antibiotic is responsible for PAS observed in the present study. Phage antibiotic combination treatments can reduce rate of resistance evolution to either phage or antibiotic or for both [73]. In the present study, U1G+CR+M & U1G+M along with colistin caused a significant 3.0 and 3.5 log decline respectively (Figure 6), which reaffirms the potential of phage cocktail to restrict bacterial bioburden in vivo and can be potentially evaluated for its efficacy in mammalian models.

5. Conclusions

Our study revealed thatthephage cocktails and colistin combination is effective in curtailing the growth of colistin resistant E. coli (U1007) both in vitro and in vivo. Machine learning tools predicted potential host cell receptors, among which, OmpC was validated as the host cell surface receptor for U1G by growth under different physiological conditions followed by the estimation of phage titers. As the lytic potential of U1G phage is enhanced by M phage in combination with colistin, Phage cocktail (U1G+M) and colistin has the potential to curtail colistin resistant E. coli in mammalian models.