# Prediction of strain level phage-host interactions across the Escherichia genus using only genomic information

Baptiste Gaborieau, Hugo Vaysset, Florian Tesson, Inès Charachon, Nicolas Dib, Juliette Bernier, Tanguy Dequidt, Héloïse Georjon, Olivier Clermont, Pascal Hersen, et al.

► To cite this version:

Baptiste Gaborieau, Hugo Vaysset, Florian Tesson, Inès Charachon, Nicolas Dib, et al. Prediction of strain level phage-host interactions across the Escherichia genus using only genomic information. Nature Microbiology, 2024, 9 (11), pp.2847-2861. 10.1038/s41564-024-01832-5. hal-04781608

HAL Id: hal-04781608

https://hal.science/hal-04781608v1

Submitted on 16 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Prediction of strain level phage-host interactions across the Escherichia genus using**
2 **only genomic information**

3
4 Baptiste Gaborieau<sup>1,2,3†*</sup>, Hugo Vaysset<sup>4,5,6†</sup>, Florian Tesson<sup>1,5†</sup>, Inès Charachon<sup>1</sup>, Nicolas Dib<sup>1</sup>, Juliette

5 Bernier¹, Tanguy Dequidt³, Héloïse Georjon⁵, Olivier Clermont¹, Pascal Hersen⁶, Laurent Debarbieux⁶,

6 Jean-Damien Ricard<sup>1,2</sup>, Erick Denamur<sup>1,7</sup>, Aude Bernheim*<sup>4</sup>

7 † These authors contributed equally

8 * Corresponding authors

9 ¹ Université Paris Cité, INSERM, UMR1137, IAME, Paris, France.

10 ² AP-HP, Hôpital Louis Mourier, DMU ESPRIT, Service de Médecine Intensive Réanimation, Colombes,

11 France.

12 ³ Institut Pasteur, Université Paris Cité, CNRS UMR6047, Microbiologie Intégrative et Moléculaire,

13 Bacteriophage Bacterium Host, Paris, France.

14 ⁴ AgroParisTech, Université Paris-Saclay, Paris, France.

15 ⁵ Institut Pasteur, Université Paris Cité, INSERM U1284, SEED, Molecular Diversity of Microbes lab,

16 Paris, France.

17 ⁶ Institut Curie, Université PSL, Sorbonne Université, CNRS, Laboratoire Physico Chimie Curie, UMR

18 168, Paris, France.

19 ⁷ Université Paris Cité, AP-HP, Hôpital Bichat, Laboratoire de Génétique Moléculaire, Paris, France.

20 Corresponding author email: aude.bernheim@pasteur.fr, baptiste.gaborieau@aphp.fr
21 Abstract

22 Predicting bacteriophage infection of specific bacterial strains promises advancements in phage therapy
23 and microbial ecology. Whether or not the dynamics of well-established phage-host model systems
24 generalize to the wide diversity of microbes is currently unknown. Here, we show that we could accurately
25 predict the outcomes of phage-bacteria interactions at the strain level in natural isolates from the genus
26 Escherichia using only genomic data (AUROC of 86%). We experimentally established a dataset of
27 interactions between 403 diverse Escherichia strains and 96 phages. Most interactions are explained by
28 adsorption factors as opposed to antiphage systems which play a marginal role. We trained predictive
29 algorithms and pinpoint poorly predicted interactions so as to direct future research efforts. Finally, we
30 established a pipeline to recommend tailored phage cocktails, demonstrating efficiency on 100 pathogenic
31 *E. coli* isolates. This work provides quantitative insights into phage–host specificity and supports the use
32 of predictive algorithms in phage therapy.
33 **Introduction**

34 Using bacteriophages (phages), viruses infecting bacteria, to combat bacterial infections in a clinical setting
35 (phage therapy), is a potential solution for addressing bacterial pathogens that are resistant to multiple
36 antibiotics¹,². Predicting which phages precisely target and infect bacteria offers potential to accelerate the
37 availability and efficacy of treatments by phages³.

38 Successful lytic infection of a bacterial cell by a phage is a multistep process involving different bacterial
39 and phage traits⁴, including bacterial hosts' phage receptors (e.g. surface polysaccharides, membrane
40 proteins) and its antiviral arsenal, as well as their viral counterparts: phage Receptor-Binding Proteins
41 (RBP, e.g. tail proteins) and anti-defense mechanisms (e.g. anti-CRISPR or anti-Restriction Modification
42 systems)⁵⁻⁸. Those traits can be highly diverse in natural populations, even among closely related bacteria
43 and phages⁹,¹⁰. Consequently, most phages are only capable of infecting hosts within a single bacterial
44 species or genus, and are even often strain-specific¹¹⁻¹⁴. Since most studies on phage-bacteria interactions
45 have focused on a few model couples (e.g. *E. coli* K-12 strain, and/or lambda or T phages)¹⁵⁻¹⁸, it remains
46 unknown to what extent the knowledge acquired on these model organisms could be sufficient to capture
47 the broad diversity of phage-bacteria interactions that are encountered in the wild.

48 Recent development of bioinformatic tools to characterize bacterial and phage traits from their genomes
49 now enables a finer characterization of diverse features involved in phage-bacteria specificity¹⁰,¹⁹⁻²¹.
50 Current algorithms typically determine which bacterial genus a given phage can infect but cannot make
51 predictions at deeper taxonomic levels²²⁻²⁵. Recent studies have used genomic traits to decipher phage-
52 bacteria interactions within different bacterial species²⁶⁻³⁰ paving the way to predict phage-bacteria
53 specificity from their genomes. In this work, we set out to predict phage-host specificity from genomic data
54 in a large and diverse set of bacterial strains and phages within the *Escherichia* genus, a major pathogen
55 worldwide³¹.
56 **Results**

57 **A collection of diverse *Escherichia* strains and phages**

58 To build algorithms capable of predicting the interaction between any *Escherichia* strain and any phage,
59 we first needed an experimental dataset of interactions between a collection of *Escherichia* strains and a
60 collection of phages. We also required that both collections – bacteria and phages – would be sufficiently
61 large and diverse, fully sequenced and genomically characterized. To our knowledge, such a dataset was
62 not available. We thus decided to build one *ex nihilo*.

63 **The Bertrand Picard collection: 403 bacterial strains, representative of the *Escherichia* genus diversity.**
64 First, we put together a collection of 403 natural isolates belonging to the *Escherichia* genus that we named
65 after the French medical bacteriologist Bertrand Picard who initiated this collection in the early 1980's
66 (Biography in **Supplementary Text 1**). This collection is an assemblage of 369 natural isolates, which
67 have been phenotypically well characterized in previous studies³²⁻³⁵ and of 34 strains that were used to
68 isolate the phages in this study. A focus when establishing this collection was to encompass *Escherichia*
69 diversity (**Figure 1A-E, Supplementary Text 2, Supplementary Table 1**).

70 We sought to characterize the traits that could potentially influence phage susceptibility in bacterial
71 genomes. We identified the bacterial surface structures known to be potential receptors for phages⁴,⁶
72 unveiling a wide variety of these structures: O-antigen (Figure 1C), type of outer-core of
73 Lipopolysaccharides (LPS) (Figure 1F), capsule serotypes from ATP-Binding Cassette (ABC)-dependant
74 capsule (Figure 1G) and Klebsiella capsule (Supplementary Table 1), and membrane proteins (Figure
75 1H, Supplementary Table 4). We next focused on intra-cellular defense, identifying known antiphage
76 systems (Supplementary Table 5). On average, an isolate encodes 8 defense system subfamilies, with
77 considerable heterogeneity among strains (i.e. from 1 to 16 defense systems) (Figure 1I). Some systems
78 are widespread (e.g. MazEF or RM type I), while many are present in only a few isolates (Figure 1J). As
79 a result, our genomic analysis of the Picard collection of *Escherichia* strains provides a large panel of
80 genomic traits relevant to investigate phage-bacteria interactions at the scale of the genus.

81 **The Antonina Guelin collection: 96 phages, covering a large taxonomic diversity of *Escherichia* phages.**
82 We then built a collection of 96 *Escherichia*-infecting virulent phages that we named after the French
83 physician and biologist Antonina Guelin who led the bacteriophage laboratory in the Institut Pasteur for 25
84 years (Biography in **Supplementary Text 3**). We isolated 94 of these phages from sewage water, collected
85 in different locations around the Paris area (France) over a decade (see methods). The two archetypal
86 virulent *E. coli* phages T4 and T7 were added to the 94 other phages. Importantly, the 96 phages were
87 isolated on 34 bacterial strains which are included within the *Picard* collection (**Figure 1A**). The strains
88 used to isolate and propagate phages cover 66% of *E. coli* phylogroups. All phages in the *Guelin* collection
89 belong to one of the three *Caudoviricetes* morphotypes (Myo-, Podo- and Siphovirus). This collection
90 covers 8 subfamilies and 19 phage genera (**Figure 1K**)³⁶,³⁷.

91 Altogether, we built and genomically characterized two reference collections of bacterial isolates and
92 phages specific to the *Escherichia* genus (see data availability). We hypothesized that those two collections
93 would be sufficiently large and diverse to infer generalizable rules governing phage-bacteria interactions at
94 the strain level in the *Escherichia* genus.
<figure><img src="image_2.png" /><figcaption>Figure 1. Generation of two large and diverse collections of *Escherichia* bacterial strains (*Picard collection*) and phages (*Guelin collection*).</figcaption></figure>

95

96

97
A-J The *Picard* collection of bacterial strains (n = 403) encompasses the phylogenetic diversity of the *Escherichia* genus and features a wide range of potential phage receptors and antiphage systems. A. Phylogenetic tree of bacterial isolates. Each species of the genus *Escherichia* and each of the eight phylogroups (A to G) and Clade I belonging to the *E. coli* species are represented by a warm colour from brown to orange. The three rings represent respectively the different 163 ST, the 93 O-antigen types and the 41 different H-antigen types. B. Frequency of the bacterial sequence types. C. Frequency of the bacterial O-antigens. D. Number of isolates by pathotype. ExPEC: Extraintestinal Pathogenic *E. coli*. InPEC: Intranintestinal Pathogenic *E. coli*. E. Number of isolates by ecological source. Human sources include feces, urine or blood samples. F. Number of isolates according to their lipopolysaccharide (LPS) outer core type. Some strains could not be typed for this trait (n = 36). G. Frequency of the bacterial capsular serotypes (ABC capsule type). Of the 171 isolates with an ABC-dependent capsule, we identified K-antigen for 101 of them, retrieving the most described *Escherichia* capsular serotypes. H. Number of variants of a putative phage receptor outer membrane protein. Out of the 12 outer membrane proteins previously identified as putative receptors for *Escherichia* phages⁶, a large diversity of variants (gene clusters with more than 99% genomic identity) has been identified. I. Number of defense system subfamilies per bacterial genome. J. Frequency of subfamilies of bacterial defense systems. K. Network representation of 1,122 *Escherichia* phage genomes. The *Guelin* phage collection (n = 96, colored dots) encompasses the diversity of the *Caudoviricetes* phages targeting the *Escherichia* genus. Their taxonomic distribution was analyzed in the context of 1,026 complete genomes of *Escherichia* phages (in gray, *Genbank* database). Each genome is represented by a node (point); edges (line) were calculated using the number of shared protein groups. Guelin collection's phages are colored by their genera, which are grouped by morphotype.

## An experimental dataset of phage-bacteria interactions

Following the establishment of the *Picard* and *Guelin* collections, we assessed the susceptibility of each strain to each phage. For that, we tested each possible phage-bacteria combination. The susceptibility of the strains to phages was assessed by plaque assay experiments (see methods) using three different phage concentrations corresponding to an approximate ratio between the number of virions and the number of bacterial cells of 10, 1, and 0.1 (Figure 2A). For each phage-bacterium couple, the interaction data points were aggregated into a single interaction score: the "Minimum Lytic Concentration" (MLC) score which corresponds to the lowest concentration of phage at which a lytic interaction was observed (see methods). Such encoding allows to score each interaction in a semi-quantitative way. (Figure 2B). The resulting phage-bacteria interaction matrix comprises 38,688 phage-bacteria interactions (Figure 2C). We provide all associated pictures of plaque assay experiments to allow exploration of this data set in future studies (see data availability).

Each phage was able to infect an average of 83 (±70) bacterial strains (20.6% of the *Picard* collection) (Figure 2D), while a bacterial isolate could be lysed by 20 (±15) phages (20.8% of the *Guelin* collection) on average (Figure 2E). Only 12 (3.0%) isolates were not infected by any of the phages. We observed significant discrepancies among groups of phages and groups of bacteria. While our knowledge of phage-bacteria interactions in *Escherichia* is essentially based on model bacterial strains (K-12, *Escherichia B*), less than half of the phages were able to infect *Escherichia B* strains and 35.4% *E. coli* K-12 (Supplementary Figure 1). Myoviruses have a significantly broader host range than Podo- and Siphoviruses (41.3% of lytic interactions on average vs. 11.1% and 9.1% respectively (One-way ANOVA: F = 114.7, *p* < 10⁻³, Figure 2F and Supplementary Text 4). The full matrix exhibits an overall nested pattern (Supplementary Text 5 and Supplementary Figure 2).
A

<figure><img src="image_3.png" /></figure>

B

<figure><img src="image_4.png" /></figure>

D

<figure><img src="image_5.png" /></figure>

E

<figure><img src="image_6.png" /></figure>

F

<figure><img src="image_7.png" /></figure>

C

<figure><img src="image_8.png" /></figure>
Figure 2. A matrix of 38,688 phage-bacteria interactions in the *Escherichia* genus reveals complex patterns.

A-B Schematics of phage-bacteria interaction outcome determination. A. Experiment set up. We investigated each phage-bacteria interaction combination by plaque assay at different phage concentrations (i.e. approximate phage:bacteria ratios of 10, 1, and 0.1). Each experiment was performed in triplicate. B. Schematic of the assignment of MLC scores. Each column corresponds to an example of how the MLC score was determined. Phage-bacteria interactions were encoded using the MLC that corresponds to the lowest concentration of the phage at which a lysis of the bacterial lawn was observed. 0 indicated that no lytic interaction was observed; 1 is lytic interaction at the highest phage titer (5 x 10⁸ PFU/mL); 2 is lytic interaction at middle phage concentration (5 x 10⁷ PFU/mL); interactions observed at the lowest phage concentration (5 x 10⁶ PFU/mL) are distinguished between 3 (individualized lysis plaque) and 4 (entire lysis of the bacterial lawn) (see methods). MLC: minimum lytic concentration. C. Phage-bacteria interaction matrix. Each column is a phage. Each line is a bacterial isolate. Colors correspond to phage viral genera and bacterial phylogroups. Phage and bacteria are distributed according to their phylogenetic proximity as determined in figure 1. Intensity of blues in the matrix corresponds to the MLC score. D. Percentage of phages infecting per bacterium. The graph represents the number of bacterial isolates stratified by the percentage of the 96 bacteriophages capable of infecting them (per 5% phage susceptibility). E. Percentage of isolates infected by each phage. The graph represents the number of phages stratified by the percentage of the 403 isolates they can infect. F. Breadth of the phage host range according to the phage morphotype. Each data point corresponds to a phage of the *Guelin* collection (n = 96). Phages are grouped by morphotype (Myovirus, Podovirus, Siphovirus). In each boxplot, the box is centered on the median and spans from the first to the third quartile. Whiskers span to the furthest datapoint within the 1.5 × interquartile range (One-way ANOVA, F=114.7, p-value=1.44×10⁻²⁶).

Overall, the global analysis of the interaction matrix shows that our broad collection of phages, characterized by their diversity, results in a lytic interaction in 20.3% of phage-bacteria pairs, at the strain level within the *Escherichia* genus.

Phage RBPs heavily correlate with phage host range

We then wanted to use this experimental dataset to predict phage-bacteria interactions. A straightforward approach would be to use the 38,688 phage-bacteria interactions in our dataset to train a machine learning algorithm. The input features taken by such an algorithm could a priori be anything that potentially influence bacterial susceptibility to phages, including adsorption factors or antiphage defense systems. However, the diversity of both adsorption factors (e.g. > 90 serotypes of O-antigen) and defense system families in the collection (>100 families) makes the number of candidate traits provided as input features considerable. This could cause machine learning to overfit their training dataset and loose generalization capability. We thus seeked to narrow down the number of input features to our models by better understanding which genomic traits govern phage-bacteria interactions in our dataset and quantifying their explanatory power.

First, we examined phage infection patterns in the matrix. We found that phages from the same viral genus (>70% nucleotide identity of the full genome length) share more similarity in their infection patterns than what would be expected for two phages taken at random (PERMANOVA, adjusted R²=0.13, *p* < 10⁻⁵)³⁸.
We further observed that the major explanatory variable of phage host range is its isolation strain (PERMANOVA, adjusted R²=0.53, *p* < 10⁻⁵). The interaction between both variables accounted for more than 60% of the variance in the distance matrix (Permutational ANOVA, adjusted R²=0.65, *p* < 10⁻⁵) (Figure 3A)³⁸. We then investigated which specific genomic traits could explain phage infection.
Comparative genomics between phages of the same genus showed that most of the core genes were highly conserved in terms of amino acid identity within each genus, while genes encoding for tailspike or tail fiber proteins were systematically hotspots of variability³⁹. Additionally, we observed that those tailspike and tail fiber-encoding genes were more conserved among phages isolated on the same strain, compared to phages isolated on different isolation strains (Figure 3B and Supplementary Figure 3). Tail fiber and tailspike genes are well known to encode for RBP⁴,⁴⁰,⁴¹. We thus hypothesized that RBP could be the major phage traits explaining the infection patterns observed in the interaction matrix, corroborated by some examples in our dataset (Supplementary Text 6, Supplementary Figure 4). The systematic identification of RBP in the 96 phages showed that the number of lytic interactions performed by a phage was not linked with the number of RBPs it encodes, but rather likely to be influenced by the targeting specificity of its RBPs (Supplementary Text 7).

<figure><img src="image_9.png" /><figcaption>Figure 3. Phage isolation strain and Receptor-Binding Proteins are major determinants of the phage host range.</figcaption></figure>

A. Phage genus and isolation strain associate with phage host range. Left: Dimensionality reduction of the phage host range using the Uniform Manifold Approximation and Projection (UMAP) algorithm⁴³. The 403-dimensional vector encoding the binary interactions of each phage with all the bacteria is reduced to a 2-dimensional vector. The color of each phage indicates its genus and the letter indicates the bacterial strain on which it was isolated and propagated. Right: Host range of the phages located within the two zoom boxes present on the UMAP graph. Phages are clustered according to their host range similarities. Shades of blue indicate the MLC score assigned to the interaction between each phage and each bacterium (Ø: No infection). B. Phage RBP are hotspots of genomic variability and are more conserved among phages isolated on the same host. Top Left: Genomic map of six representative phages of the *Vectrevirus* genus generated using Clinker⁴⁴. The six phages were isolated on four different strains indicated by a letter. Genes are colored according to their functional annotation: either RBP (blue) or other (gray). The stroke between each pair of genes at the same locus indicates the percentage of aminoacid sequence similarity. Bottom Left: Position-wise entropy of the genomic alignment. A multiple whole genome alignment is performed including the six phage genomes and the Shannon entropy is computed at each position of the alignment based on the distribution of the
nucleotides. The lower the entropy, the stronger is the conservation of nucleotides at a given position in the multiple sequence alignment. Right: Host range of the six phages presented in the genomic map.

Overall, our general statistical analysis supports the idea that phages isolated on the same host exhibit similar infection patterns because they tend to share similar RBP, as previously suggested in *Salmonella*<sup>42</sup>. As such, the diversity of isolation strains plays a fundamental role when selecting phages with a diverse set of RBPs. This implies that the choice of the isolation strains is an important parameter to take into account when isolating phages, especially in the context of phage therapy.

## 217 **Adsorption factors drive bacterial susceptibility to phages**

We investigated which bacterial traits exhibit the strongest association with the bacterial infection patterns, notably trying to evaluate the relative importance of adsorption and intracellular defense. To do so, we assessed the association of bacterial adsorption factors and defense systems with the probability of being infected by each phage. We fitted one Generalized Linear Mixed-Model (see method) (GLMM) per phage, aiming at explaining the lytic interactions of its corresponding phage based on both the adsorption factors and the defense systems subfamilies<sup>14,45</sup>. A total of 30 bacterial adsorption factors significantly associated with phage lytic interactions whereas only 2 defense systems did the same (**Figure 4A**). We also followed a complementary statistical approach which gave similar results (**Supplementary Text 8**). This suggests that bacterial adsorption factors were the best variable to explain the overall infection patterns, mirroring what is observed in phages.

We then sought out to characterize more precisely which adsorption factors would associate the most with phage infection in the matrix. To do so, we fitted a distinct GLMM on each individual phage infection pattern, taking only the bacterial adsorption factors as covariates. Here, we identified a total of 27 adsorption factors significantly associated with the lytic interactions of 24 phages in the *Guelin* collection. The most frequent adsorption factors yielded by the linear models were mostly surface polysaccharide-related traits (K-serotype, LPS outer core or O-antigen serotype). Outer membrane proteins – frequently identified as a phage receptor in several studies – were not significantly associated with phage infection in our dataset<sup>6,46</sup> (**Figure 4B**). Importantly, specific adsorption factors such as the *Klebsiella* capsule serotype exhibited a significant but weak quantitative association with the infection patterns in the matrix (PERMANOVA, adjusted R<sup>2</sup>=0.02, *p*<10<sup>-3</sup>). In spite of their weak quantitative association with infection patterns, those traits are worth investigating as they highlight local, non trivial, patterns in the matrix (e.g. **Supplementary Figure 5**).

## 240 **Antiphage systems play marginal role in predicting interactions**

As we were surprised by the limited explanatory power of antiphage systems on the host range in our analyses, we further investigated their impact on phage-bacteria interactions. First, we evaluated if bacteria encoding more defense systems would be protected against more phages. We found no evidence of such association (Linear Mixed Model (LMM), $\beta_{Num\,DS}$= -0.22, 95% credible interval = [-0.74, 0.34], pMCMC=0.41) (**Figure 4C**)<sup>45</sup>. We further hypothesized that bacterial antiphage systems function once a phage has adsorbed and is starting the replication cycle. This would imply that an increase in the number of defense systems encoded by a strain would not impact the number of infecting phages, but rather would make infecting phages less virulent (i.e. with a decreased MLC score). We observed a low but significant correlation between the number of defense systems encoded in a strain and its average MLC upon phage
infection (taking only the lytic interactions into account) (LMM, $\beta_{Num\,DS}$ = -0.04, 95% credible interval = [-0.05, -0.01], pMCMC=4×10⁻³) (Figure 4D, Supplementary Text 9, Supplementary Figure 6). Those results strongly suggest that the number of defense systems of a bacterial strain does not influence the probability of infection by phages in our interaction dataset but could slightly impact the level of susceptibility to already infecting phages. Overall, these results show that bacterial defense systems are not required to predict the phage-bacteria interactions for the *Escherichia* genus and can be removed from the set of candidate traits provided as input features to our models.

<figure><img src="image_10.png" /><figcaption>**Figure 4.** Bacterial adsorption factors are the major determinants of phage-bacteria interactions in contrast with defense systems which marginally reduce virulence of infecting phages.</figcaption></figure>

A. Number of bacterial adsorption factors and defense systems significantly associated with phage lytic interaction. A distinct binomial GLMM is fitted on the lytic interactions performed by each phage. GLMM are provided with the bacterial adsorption factors and the defense systems subfamilies as covariates. We then count how many times each bacterial trait (adsorption factor or defense system) significantly associates with lytic interactions by each phage. We consider that the association is significant if pMCMC is lower than 10⁻³ (see methods)⁴⁵. B. Number of times each individual bacterial adsorption factor associates significantly with phage lytic interaction. As in A, one binomial GLMM is fit on the interactions of each individual phage with the bacterial phylogeny as a random effect⁴⁵. Each model takes only the bacterial adsorption factors as covariates. K-type: ABC dependent K-antigen serotype. LPS OC: LPS Outer Core. O-type: O-antigen serotype. C. Number of infections undergone by each bacterial isolate as a function of the number of detected antiviral defense systems. Each data point corresponds to a bacterial strain in the *Picard* collection (n = 403). The best regression line obtained by fitting a LMM using the number of defense systems as a covariate and taking the core bacterial phylogeny into account as a random effect is shown⁴⁵. $\beta(Num\,DS)$: mean of the posterior distribution of the coefficient associated with the number of defense systems. D. Average MLC observed for each strain as a function of the number of detected antiviral defense systems. As in C, each data point corresponds to a bacterial strain in the *Picard* collection (n = 403) and the best fit line obtained by fitting a LMM taking the number of defense systems as a covariate and the core bacterial phylogeny as a random effect is represented. Only lytic (MLC > 0) interactions are included when fitting the LMM.

## 277 Adsorption factors are sufficient to predict interactions

Identifying the bacterial and phage genomic traits allowed us to dramatically narrow down the number of traits that we should provide as input features to a model aiming at predicting phage-bacteria interactions at the strain level in the *Escherichia* genus. Given the results of the statistical inference, we hypothesized that the adsorption factors alone would be sufficient to predict the phage-bacteria interactions. To test this hypothesis, we trained machine learning models to predict the phage-bacteria interactions as a binary classification task (i.e. lytic vs. non lytic interaction) (see methods). Phage-specific models were evaluated on their capacity to generalize to previously unseen bacterial strains and achieved an Area Under the ROC
curve (AUROC) of 85% on average (Figure 5A). When considering predictions for each of the 38,688 phage-bacteria interactions as a whole, an overall AUROC of 86% (Average Precision = 60%, F1 score = 60%) was obtained on Group 10-Fold Cross-Validation (Figure 5B & D). This remains true for bacteria with as few as 10% or less of lytic interactions which indicates good generalization ability even in the case of class imbalance (Figure 5C)⁴⁷. The prediction performance we achieve shows that the adsorption factors are sufficient to accurately predict phage-bacteria interactions in our dataset. Our detailed understanding of which phage-bacteria interactions are correctly predicted provide a roadmap for further research as it will allow us to focus on the 3,379 false negatives and 2,922 false positives that are not captured by our current genomics-based predictive model over a total of 38,688 interactions (Figure 5D).

<figure><img src="image_11.png" /><figcaption>Figure 5. Adsorption factors are sufficient to accurately predict phage-bacteria interactions, even for phages with very few lytic interactions.</figcaption></figure>

A. Distribution of the prediction performance among phage-specific models. One classification model is fitted for each individual phage and evaluated on Group 10-Fold Cross-Validation. The performance on the validation sets of each phage-specific model is measured by the AUROC metric and averaged across all cross-validation folds. One AUROC score is computed per bacterial strain to measure the ability of each model to generalize to previously unseen bacteria. B. Overall prediction performance of the predictive algorithm across the 96 phages. The predictions of each phage-specific model on the validation set at each cross-validation fold are aggregated together into one set of 38,764 predicted interactions. Prediction performance of the obtained model is assessed using the ROC curve (blue) and the Precision-Recall curve (green). C. Distribution of the prediction performance among phage-specific models as a function of the phage number of lytic interactions. The performance of each phage-specific model (n =
96) on the validation sets of each cross-validation fold is averaged. It is plotted against the number of lytic interactions performed by its corresponding phage, which gives the amount of class imbalance in the classification problem. D. Prediction matrix obtained using adsorption factors and core phylogeny only. One classification model per bacteriophage is trained to predict the phage-bacteria interactions bacteriophage from the interaction matrix. Each model takes as input the same bacterial adsorption factors as well as the core phylogeny of the bacteria. They are trained and evaluated on Group 10-Fold Cross-Validation. Predictions obtained on the validation set at each round of the cross-validation procedure are compared with the ground truth and are reported in this prediction matrix. Blue: True Positive. Red: False Positive. Gray: False Negative. White: True Positive.

## 314 Recommending phage cocktails to target pathogenic strains

We next evaluated whether our predictive models for phage-bacteria interactions could be used in a framework of phage therapy where a combination of several virulent phages is assembled into a phage cocktail to target a specific pathogenic bacterial strain (hereafter referred to as the “query” strain)⁴⁸. We aimed at designing a machine learning-based phage cocktail recommender system to target pathogenic *E. coli* strains based on bacterial genomic traits only. We then wanted to evaluate the performance of this algorithm *in vitro* on a test collection of previously unseen *E. coli* pathogenic isolates (Figure 6A).

We decided to use a cocktail size of three phages offering a good trade off between efficiency and parsimony (Supplementary Text 10, Supplementary Figure 7). We then went on to design a recommender system that goes sequentially through four steps and stops whenever three phages are recommended (Figure 6A, Supplementary Figure 8). The four steps are detailed in the methods section. The rationale behind this pipeline organization of the recommender algorithm was to prioritize steps that were empirically measured as more precise but often less exhaustive over steps that are more exhaustive but less precise.

We then assessed the capacity of our recommender system to generalize to a new collection of 100 previously unseen *E. coli* isolates responsible for Ventilator-Associated Pneumonia (VAP) (Figure 6A-B, Supplementary Table 6). Phage cocktails were recommended based solely on the genomic traits of these 100 strains. A total of 18 different phage cocktails were recommended for the whole test collection (Supplementary Table 7). We observed that the recommended cocktails could be gathered into two categories based on the steps of the pipelines at which phages were recommended. First, some cocktails (n=15, recommended to 24 test bacteria) contain phages recommended at early stages of the pipeline (denoted as “tailored” in the rest of the text) while the rest of the cocktails (n=3, recommended to 76 bacteria) are recommended at late stages of the pipeline (denoted as “generic” in the rest of the text) (Figure 6C). This behavior is a direct consequence of the pipeline organization of our recommender algorithm and reflects the fact that some bacteria – those that recommended tailored cocktails – are easily handled by the pipeline whereas others are handled with more difficulty and tend to recommend generic cocktails. As a consequence, tailored cocktails sometimes contained highly specific phages (e.g. T145_P2 which infects 13.1% of the original collection) often recommended early in the pipeline compared with generic cocktails which only contained broad range phages often recommended downstream in the pipeline (Supplementary Figure 9). To evaluate the performance of our algorithm, we compared the performance of the tailored and generic cocktails with an uninformed baseline cocktail. The latter corresponds to the most covering triplet of phages, i.e. the cocktail of the three phages which together lyse the largest number of bacteria in the original interaction matrix. This baseline cocktail yielded productive lysis on 63% of the original Picard collection.
We then evaluated experimentally the tailored, generic and baseline cocktails by plaque assay experiments. No antagonistic or synergistic phage-phage interaction was observed in this experimental setup (Supplementary Text 11). To evaluate the performance of the cocktails, we first compared the overall percentage of lysed bacteria for each method. The baseline and the generic recommendations had similar percentage of success (81.00% and 78% lysed bacteria) whereas the tailored cocktails had a higher success rate on the bacteria it was recommended to (91.67% lysed bacteria) (Figure 6D). Second, we evaluated the MLC of the full cocktail. The tailored cocktails significantly outperformed both the baseline (two-sided Mann-Whitney test: U=399.5, *p*=0.02) and the generic cocktails (two-sided Mann-Whitney test: U=1171.5, *p*=0.02, Figure 6E). Overall, our algorithm successfully recommended tailored phage cocktails based only on genomic data. The recommended cocktails could complement a broad baseline approach in a typical clinical use.

<figure><img src="image_12.png" /><figcaption>Figure 6. A recommender system allows the design of phage cocktails targeting previously unseen pathogenic bacteria based on their genomic traits.</figcaption></figure>

A. Schema of the experiment. The recommender algorithm trained on the original interaction matrix is used to predict the interaction between a query strain (i.e. from a “test” collection of 100 unseen *E. coli* pathogenic isolates) with any phage of the *Guelin* collection based solely on its genomic traits. A cocktail of three phages is recommended to the query strain. The recommended cocktail is then evaluated by challenging each strain from the test collection with its recommended cocktail as well as with a baseline cocktail (three most covering triplet of phages) at six different phage concentrations (from 5 x 10³ to 5 x 10⁸ PFU/mL). B. Phylogenetic tree of the strains. Core phylogeny of the 100 test strains (black triangles) is computed jointly with the 403 strains of the *Picard* collection. C. Number of *E. coli* strains in the test collection to which each of the 18 different cocktails was recommended. Recommended cocktails (named with a letter from A to R, *n* = 18 distinct cocktails) are categorized into either «generic» or «tailored» according to the steps of the pipeline at which they were recommended. D. *In vitro* success rate of the recommended cocktails on the test collection. Success rate is the percentage of lytic interactions performed by each category of cocktail over the total number of times such a cocktail category was recommended. Cocktails are grouped according to their category (Baseline, Generic, Tailored). Each cocktail (*n* = 100), is tested on three biological replicates. E. Binned MLC
obtained *in vitro* on the test collection. The MLC score (ranging from 0 to 6) is binned into four ordered categories: Ø («No infection», MLC=0), L («Low», $1 \le \text{MLC} \le 2$), M («Medium», $3 \le \text{MLC} \le 4$) and H («High», $5 \le \text{MLC} \le 6$). Each cocktail (n = 100) is tested on three biological replicates. Asterisks indicate a statistically significant difference among cocktail groups (two-sided Mann-Whitney test, *p*(Tailored vs. Baseline) = 0.019, *p*(Tailored vs. Generic) = 0.034 and *p*(Generic vs. Baseline) = 1.0), n.s.: not significant. MLC: minimum lytic concentration.

## 380 Discussion

Our work demonstrates that it is feasible to predict phage-host specificity from genomes with good accuracy. Among 38,688 interactions between 403 *E. coli* isolates and 96 phages, we understand 65% of the variance in phage host range based on the identified phage genomic traits. We predict phage-bacteria lytic interactions from bacterial genomic traits with an AUROC of 86%. We further show that such predictions can be harnessed to design tailored and efficient phage cocktails.

Our predictions could be improved and refined (Supplementary Text 12).
Despite limitations, our findings provide a clear picture on the genomic factors governing phage-host specificity in our collection of natural *E. coli* isolates, where adsorption factors play a major role, while antiphage systems have marginal impact. Unlike previous studies on *Vibrio crassostreae*¹¹,¹², our data doesn't show that antiphage systems shape the phage host range. This may stem from the broader evolutionary scope of our collection. The broader evolutionary scale allows a higher diversity of adsorption factors while the defense turnover is already high in smaller phylogenetic distances⁴⁹. We showed that more defense systems are associated with lower virulence of infection but the effect size is low. This could be explained by the "pan-immune system" hypothesis (natural bacterial communities are heterogeneous and individual strains encode only a small subset of all the defense potential of the community). This could also suggest alternative ecological roles for antiphage systems, such as being involved in competition between mobile genetic elements⁵⁰,⁵¹. Whether this holds true for other bacterial genera, remains to be determined.

Our findings not only quantify key drivers of phage-bacteria interactions but also provide a roadmap for further understanding. Both the depth and breadth of the matrix and the genomic and phenotypic characterization allow to map the existing knowledge efficiently and can thus orient future studies to the interactions not explained by our current models. This systematic approach could be applied to diverse bacterial species, including other ESKAPE pathogens to advance our understanding.

Finally, we used predictions to inform phage therapy applications by demonstrating the need to isolate phages on diverse strains and by training a phage cocktail recommender system. Our pipeline, along with other recent works²⁹, allows clinicians to input a pathogenic strain's genome, receiving tailored cocktail recommendations for heightened success rate and a stronger lysis. This algorithm, leveraging genomic data, streamlines phage selection, enhancing therapy efficiency within a personalized medicine framework.
408 Methods

409 Culture media and solutions

410 Strains were cultivated at 37°C on LB Lennox medium (agar or broth) (Becton Dickinson, USA, ref. 240110
411 & 240230). The thickness of the agar plates was standardized using an automatic pump dispensing 50 mL
412 per 120mm x 120mm petri dish on a flat surface. Phosphate-buffered saline (PBS) 1x was prepared by
413 dissolving commercial tablets (Sigma-Aldrich, USA, ref. P4417) in 200 mL of deionized water and
414 sterilized by autoclaving. All the batch media were prepared within a few days to ensure the most
415 homogeneous condition.

416 Picard collection (bacterial strains)

417 The *Picard* collection includes 369 strains, previously gathered to be representative to the phylogenetic
diversity of the *Escherichia* genus and described by Galardini et al.<sup>34,35</sup>, and the 34 strains used to isolate
418 and propagate the phages were also added to this collection. All the details of the 403 strains are available
419 in the **Supplementary Table 1**. All assembled genomes are available on Figshare under the DOI
420 10.6084/m9 figshare.25941691.
421

422 **ColoColi** collection (bacterial strains)

423 To test the validity of our predictive phage cocktail algorithm, in a context compatible with phage therapy
424 settings, we used a collection of 100 pathogenic strains of *E. coli*. These strains originated from the *ColoColi*
425 collection, which have been previously described⁵². In brief, *ColoColi* gathered 310 strains mainly from
426 Ventilator Acquired Pneumonia (VAP) from 14 different hospital centers in France. It includes in majority
427 pathogenic strains (*n* = 263) and some colonizing strains (*n* = 47). Half of the strains used to isolate and
428 propagate phages of the *Guelin* collection came from the *ColoColi* collection. The test strains set were
429 chosen to be as diverse as possible in terms of phylogeny and geographical origin. Sampling of this
430 collection led to the selection of 100 strains, excluding the phage host strains, and representing 7
431 phylogroups of *E. coli*. No phylogroup H strains or cryptic clades were present in the *ColoColi* collection.
432 These strains represented 60 different STs (26 of which were represented in the *Picard* collection), had 32
433 different O antigens and 24 H antigens. STs that were typical of VAP such as ST127 or ST141 were kept
434 in the test strains. All the details of the 100 strains are available in the **Supplementary Table 6**.

435 Guelin collection (bacteriophages)

436 Phages were isolated from samples of sewage water collected over 10 years from several wastewater
437 treatment plants of the urban region of Paris (France). Phage isolation was previously performed on a
438 specific host with an enrichment step as described by Debarbieux et al.⁵³. Isolation host strains were all *E. coli* isolates, including clinical isolates from VAP (*n* = 21), urinary tract infection (*n* = 2), and intestinal
439 strains from inflammatory bowel disease patients (*n* = 7), as well as EAEC (*n* = 1), commensal strain of
440 mice (*n* = 1) and lab strains (K-12 and *E. coli* B, *n* =2). All phage lysates were filtered (0.2 µm) before use.
441 By convention, phage titer was determined on the host strain used to isolate and propagate the phage. Phages
442 closely related to other ones were deliberately excluded if they shared a nucleotide identity upon 98% unless
443 they displayed obvious phenotypic differences of their host-range. Ultimately, out of 128 isolated phages,
444 among which few were previously characterized⁵⁴⁻⁶⁰, 96 unique phages were obtained. All the details of the
445 96 bacteriophages are available in the **Supplementary Table 2**.
446

447 Assembling phages in cocktails
Three phages which have been predicted by the recommender algorithm to be able to lyse the query strain were assembled in a cocktail, and compared to a reference "Baseline" cocktail, composed of the three phages whose combined cover the wider host range of the *Picard* collection (LF82_P8, CLB_P2, LF73_P1). For both baseline cocktail and predicted cocktails, the phages have been tested individually (phage 1, phage 2, phage 3), in two-by-two combinations (phages 1 & 2, phages 1 & 3, phages 2 & 3) or all-together (phages 1, 2 & 3) following the same plaque assay experimental method described below. Each combination of two or three phages was a mix of an equivalent quantity (Plaques Forming Unit, PFU) of each phage, and so as to obtain phage solutions with final theoretical concentrations of 5 x 10⁸ PFU/mL. The quantity for each phage was calculated from the concentration of the phage determined on its own isolated strain. The composition of phage cocktails was detailed in the **Supplementary Table 7**.

### Evaluating phage-bacteria interaction outcomes: plaque assays experiments

The infectivity of a given phage on a given bacterial strain was quantified by plaques assays. The phage-bacteria interaction was assessed by adding a 2 µL drop of phage solution on a bacterial lawn on an agar plate. Each bacterial strain was overlaid on agar plates, from a fresh liquid medium culture, that reached the exponential growth phase (OD<sub>600nm</sub> 0.25 ≈ 5 x 10<sup>7</sup> CFU/mL). All the phage solutions titers were normalized at 5 x 10<sup>8</sup> PFU/mL in PBS. Each phage against each bacteria was tested at different phage concentrations. The multiplexing and standardization of phage deposit was performed by automation with a Viaflo 96 device (Integra Bioscience, Swiss). Agar plates were incubated overnight at 37°C. The reading of the results was standardized. All plates were digitized with a high-resolution scanner after 16 hours of incubation and reviewed manually using a custom interface built on top of the Napari interactive viewer (v0.4.16)<sup>61</sup>. The result of each interaction was encoded according to defection or not in the bacterial lawn. We encoded each phage-bacteria interaction using the "Minimum Lytic Concentration" (MLC) score, which corresponds to the lowest concentration of the phage at which a lytic interaction is observed. We consider that an interaction is lytic whenever we observe individual lysis plaques or full clearing at any phage concentration. Individual lysis plaques ascertain the lysis of the bacterial population with production of new virions. Clearing of the bacterial lawn at high phage concentration could result from productive lysis of the bacterial population by the phage, or from another mechanism such as lysis from without, or abortive infection. The average MLC upon infection was computed by taking only the lytic interactions of a given bacterial strain and averaging the MLC score across all infecting phages.

Three independent replicates were performed to assay the interactions of the 38,688 phage-bacteria couples (403 isolates from the *Picard* collection versus the 96 bacteriophages from the *Guelin* collection), at three different phage concentrations of 5 x 10<sup>8</sup> PFU/mL (replicates R1, R2 and R3), 5 x 10<sup>7</sup> PFU/mL (replicates R2 and R3), 5 x 10<sup>6</sup> PFU/mL (replicates R1, R2 and R3) and 5 x 10<sup>4</sup> PFU/mL (replicate R1). The MLC score is null in case of non-lytic interaction, and ranges from 1 (lytic interaction at the highest phage titer) to 4 (uncountable number of lysis plaques at 5 x 10<sup>6</sup> PFU/mL). Such encoding allows to score each interaction in a semi-quantitative way. The final score represents the average MLC of the three independent replicates. The outcome of interaction at 5 x 10<sup>4</sup> PFU/mL was not taken into account in the calculation of the MLC score because it was not verified by a replicate. We measured that 98.3% of the interactions – regarding any phage/bacteria/phage concentration triplet – yielded consistent annotated scores across triplicates (either same score assigned to the three replicates (88.6%) or same score assigned to 2/3 replicates (9.7%)).

Three independent replicates were performed to assay the interactions of the phages of the cocktails (alone, or in combination) with the 100 isolates from the *ColoColi* collection, at six different phage concentrations
491 of 5 x 10⁸ PFU/mL, 5 x 10⁷ PFU/mL, 5 x 10⁶ PFU/mL, 5 x 10⁵ PFU/mL, 5 x 10⁴ PFU/mL and 5 x 10³ PFU/mL. The MLC score is null in case of non-lytic interaction, and ranges from 1 (lytic interaction at the highest phage titer) to 6 (uncountable number of lysis plaques at the lowest phage titer). All the interaction results are available on GitHub (https://github.com/mdlmparis/coli_phage_interactions_2023) and on the VHRDB⁶² (https://viralhostrangedb.pasteur.cloud)
492
493
494
495

496 **Genome sequencing**

497 Phage nucleic acid was extracted by a conventional molecular method. Briefly, after a pre-treatment of the lysate with DNase (Sigma-Aldrich, USA, ref. D4527-20KU) and RNase (Thermo Scientific, USA, ref. EN5031), capsid was disrupted by adding proteinase K (Eurobio Scientific, France, ref. GEXPRK01-B5) and SDS (Sigma-Aldrich, USA, ref. 05030), and the nucleic acid retrieved by phenol-chloroform-acid-isoamyl alcohol 25:24:1 (Sigma-Aldrich, USA, ref. 77617) separation followed by a sodium acetate-ethanol (Sigma-Aldrich, USA, ref. S7899) precipitation step. Bacterial nucleic acid was extracted using the EZ1 DNA tissue kit (Qiagen, Germany, ref. 953034) on EZ1 Advanced XL (Qiagen, Germany) according to the manufacturer's recommendations. The libraries were constructed using the Nextera DNA prep Kit (Illumina, USA, ref. 20018705) according to the manufacturer's protocol. A pair-end sequencing (300 bp) was performed on a MiniSeq platform (Illumina, USA).
497
498
499
500
501
502
503
504
505
506

507 **Genome analysis**

508 Read quality was checked using FastQC⁶³. Adapters and quality reads were trimmed by TrimGalore. Reads 509 were *de novo* assembled with SPAdes v3.15.4⁶⁴ in one single contig. Ambiguities were resolved for phage 510 genomes with the Bandage visualizing tool. We used Prokka v1.14.5 (--kingdom Viruses)⁶⁵ to annotate the 511 genomes. All the genome sequences have been deposited on the ENA under the Bioproject ID: 512 PRJEB69283.
512

513 **Phylogenetic analysis**

514 *In silico* typing of the bacterial isolates was carried out using SRST2 for the sequence type (ST) according 515 to the MultiLocus Sequence Typing (MLST) Warwick⁶⁶. Phylogroup was confirmed by the Clermont 516 Typing software⁶⁷. Bacterial genomes were annotated with PanACoTA v.1.3.1⁶⁸. The phylogenetic tree was 517 built on the core genome of the strains using the PanACoTA v1.3.1⁶⁸ pipeline (pangenome, coregenome, 518 align and tree) using the fasttree option. The trees were annotated with iTOL v5⁶⁹ online tool. 519 Phage were sorted by family and genus based on whole genome BLASTn using Viridic v1r3.6⁷⁰ compared 520 against the latest release of the ICTV classification VMR_MSL38_v2 (September 2023). Virus assignment 521 into genera (≥70% nucleotide identity of the full genome length) and species (≥95% similarities) ranks 522 follows the International Committee on Taxonomy of Viruses (ICTV) genome identity thresholds. 523 For the network analysis, distance between phage were calculated by vConTACT2 v0.11.3, using both viral 524 genomes of this study as well as the 1,026 *Escherichia* phage complete genomes from *Genbank* (Organism 525 name containing “Escherichia” OR “Coli” in all viruses genomes taxid 10239). The phage taxonomic 526 network was then visualized and annotated with Cytoscape software v3.9.1⁷¹.
526

527 **Detection of genomic traits from the genomes of phages and bacteria**

528 O-antigen and H-antigen where determined using ECTyper v1.0.0⁷². Antiphage defense systems were 529 detected using DefenseFinder v1.2.0 using DefenseFinder models v1.2.3 with default parameters¹⁰, 530 CapsuleFinder v1.1 with default model for the ABC-dependent capsule except for the parameter 531 inter_gene_max_space which was set to 20⁷³. We also detected and typed the ABC-dependent capsules to
known ones using custom MacSyFinder models (available in the github repository of the project)¹⁹,²⁰. Those models were created by clustering (at 80% identity and 80% coverage) all the genes present between the boundaries of ABC capsules detected by CapsuleFinder in 1843 *Escherichia coli* complete genomes from RefSeq using mmseqs2 v13.4511¹⁷. Colocalization of the different clusters was then used to form capsule type definitions and the sequences were used to create custom HMM profiles. HMM profiles and MacsyFinder models are available (see data availability). The attribution of known capsules type was done using those models on reference K-antigen serotype loci. The *Klebsiella* capsular serotypes were determined using Kaptive v.2.0.6 with default parameters and the Kaptive reference database for *Klebsiella* K loci⁷⁵. Phage receptor binding proteins were detected using PhageRDBDetect²¹ v2.1.3 and manual curation.

## Methods regarding Nestedness and Modularity of the interaction matrix

The nestedness temperature (T) and the modularity score (Q) of the interaction matrix were assessed using respectively the BINMATNEST algorithm⁷⁶ and the BRIM algorithm⁷⁷. Null distributions for T and Q were obtained by generating n=1,000 random matrices with identical shape and average connectance (probability that a given point in the matrix yields a lytic interaction) as the original interaction matrix. T and Q scores were assessed for each random matrix and the p-value for T and Q in the original matrix was computed from the null distribution.

## Statistical inference regarding phage genomic traits

To investigate differences in host range breadth among the three phage morphotypes present in the *Guelin* collection, we performed an ANOVA followed by Tukey honest significance test.
To explore the links between phage traits and the phage host range in a multivariate approach (i.e. taking all phage-bacteria interactions into account at once), we built a distance matrix for each phage-phage pair based on their respective binary infection vector (403 components for each phage). Since most of the phage-bacteria interactions are non-lytic, the Jaccard distance was used as a metric in order to only take lytic interactions into account when comparing two phages. The obtained pairwise distance matrix was then used to (i) perform exploratory data analysis using dimensionality reduction with the UMAP algorithm⁴³ and (ii) partition the distance matrix between diverse sources of variations using PERMANOVA³⁸. In brief, PERMANOVA is a nonparametric alternative to ANOVA that allows to partition variance in a distance matrix between diverse groups defined by a covariate of interest (e.g. the phage genus), and to test if the fitted coefficients are significantly different from zero using permutations of the elements in the distance matrix. We performed PERMANOVA as implemented in the *adonis2* function of the vegan R package (v2.6-4)³⁸.

To identify which genes are responsible for encoding RBPs, we first generated a set of candidate genes using the algorithm published by Boeckaerts et al.²¹. We then manually curated this list of candidate genes using several approaches. For phage subfamilies that contained a well-studied phage (e.g. *Tevenvirinae* and T4), we exploited the conservation of gene synteny. For the rest of the phages, we performed (i) domain annotation using HMMER⁷⁸ (version 3.3.2) and the PFAM-A⁷⁹ database (version 34.0) and (ii) structure prediction using ESMFold⁸⁰ to identify if a given protein was likely to function as a RBP. HMMER hmmsearch function was run with the --cut_ga option. Functional annotation and structure prediction also allowed us to determine whether a given RBP was a tailspike (typically containing a beta-helix domain) or a tail fiber (long, unstructured monomer)⁴,⁴⁰. The accurate identification of the RBP encode in each phage genome allowed us to investigate host range similarity between three phages, T145_P2, AL505_P3 and MT1B1_P3. RBP structural similarity analysis was performed using the Foldseek Search server⁸¹. Precision
575 (TP / (TP + FP)) and Recall (TP / (TP + FN)) were used to evaluate the quality of the predictions of the
576 infection patterns using the “Bacterial strains are infected by T145_P2, AL505_P3 and MT1B1_P3 if and
577 only if they encode a K1 ABC-dependent capsule”.

578 Statistical inference regarding bacterial genomic traits

579 The relative importance between bacterial adsorption factors (surface polysaccharides, outer membrane
580 proteins) and intracellular defense systems to explain phage infection was assessed using a univariate
581 approach (i.e. by analyzing the infection patterns of one specific phage at a time). We chose to use Linear
582 Mixed-Effects models to perform this analysis. LMM allows to perform inference in cases where
583 observations are not independent, e.g. when they come from a split-plot design or when they share an
584 evolutionary history as is the case in this study. To not take this phylogenetic non-independence between
585 bacteria could lead to inflated type I error rate (more false positive detections of significant associations).
586 As the response variable is binary (non lytic vs. lytic interaction), we used a Generalized Linear-Mixed
587 Model (GLMM) which extends the framework of the linear mixed-effects model to non normally
588 distributed response variables. We fitted one binomial GLMM per phage using the MCMCglmm R package
589 v2.35<sup>45</sup>. To do so, we included both the bacterial adsorption factors and the defense systems as covariates.
590 Only traits with more than 8 observations were included in the models in order to avoid model
591 overparameterization. The inverse relatedness matrix computed from the core phylogenetic distance matrix
592 was also included in each model as a random effect to account for phylogenetic non-independence among
593 bacteria. Each phage specific model was fitted for a total of 500,000 iterations, a burning time of 50,000
594 iterations and a thinning value of 250 iterations. Again, each model was fitted using an inverse Gamma
595 prior (V=1 and v=2×10<sup>-3</sup>) and convergence was assessed by checking the effective sample size and
596 parameter autocorrelation. The significance of a statistical association was assessed using the provided
597 pMCMC associated with each trait. pMCMC is the probability of the coefficient (e.g. the slope of the
598 regression line) being 0 under its estimated posterior distribution. We chose a significance threshold of 10<sup>-3</sup>.
599

600 We also followed a complementary approach which relied on multivariate statistical models. The rationale
601 behind this approach was to treat infection patterns in the matrix as whole, instead of separating phages
602 during the analysis as was in the univariate approach. We fitted 3 PERMANOVA models: the first one took
603 both adsorption factors and defense systems as covariates (“full” model), the second took only adsorption
604 factors as covariates and the third one took only antiphage defense system subfamilies as covariates. All
605 three models using the vegan R package (v2.6-4) were fitted for 10,000 iterations and adjusted R² was
606 estimated using the sequential sum of squares. Phylogenetic non-independence of bacteria was taken into
607 account by including in the models eight covariates which were computed from the core phylogenetic tree
608 of bacteria using Principal Coordinates Analysis (PCoA). Those eight covariates were introduced as the
609 first covariates of the models so that the phylogenetic signal they encoded was taken into account before
610 the effect of additional covariates (e.g. adsorption factors or defense systems) was estimated. Finally, the
611 three models were compared based on the Akaike Information Criterion (AIC) which values both goodness
612 of fit and parsimony (number of degrees of freedom) of a statistical model. The lower the AIC, the better
613 is the model fit given a number of parameters.

614 To further investigate which adsorption factors could be explain the infection by each individual phage, we
615 fitted another GLMM per each phage taking all the bacterial adsorption factors as covariates, the bacterial
616 core phylogeny as a random effect and trying to explain its corresponding phage infection patterns. Each
617 model was fitted using an inverse Gamma prior (V=1 and v=2×10⁻³), 500,000 iterations, a burning time of
618 50,000 iterations and a thinning parameter of 250 iterations. Again, traits with less than 8 observations were
619 lumped to avoid model overparameterization and a significance threshold of 10⁻³ was chosen.

620 Finally, we estimated the statistical association between either the number of defense systems encoded in
621 each bacterial, the genome size of the strain or the number of AMR genes it encodes, and the number of
622 lytic interactions it undergoes or its average MLC upon infection. Average MLC upon infection is computed
623 by taking all the lytic interactions observed for a strain and averaging the MLC of each infecting phage.
624 This metric measures the “strength” of an infection given that a lytic interaction occurs, regardless of
625 whether the bacterial strain is infected by many phages or not. Each statistical association was estimated
626 using a simple LMM since the target variable (either the number of lytic or the average MLC upon infection)
627 is real valued and is normally distributed. Each model was provided with the inverse relatedness matrix
628 computed from the bacterial phylogeny as described in previous paragraphs and was fitted using for
629 1,000,000 iterations, with a burning time of 50,000 and a thinning parameter of 250 iterations. Model
630 convergence was checked as described earlier. Estimated posterior mean for the intercept and slope were
631 used to draw the best fit line shown in **Figures 4C, 4D** and in **Supplementary Figure 6**.

632 Training algorithms to predict phage-bacteria interactions

633 After having identified genomic traits which overall associate with the infection patterns in the matrix, we
634 sought to train predictive algorithms for the phage-bacteria interactions. We chose to predict a binary
635 response variable *i.e.* whether a given phage-bacterium couple will give a lytic interaction or not as this
636 problem formulation makes it simpler to design and evaluate predictive algorithms. For each phage, four
637 classification models were trained: two of them were Logistic Regression models with different
638 hyperparameters and the other ones were Random Forests with different hyperparameters. The rationale
639 behind this choice of models was to select one class of simple models (Logistic Regression) and one class
640 of models with higher capacity (Random Forest). More precisely, the four models were (i) Logistic
641 Regression with L2 penalty, (ii) Logistic Regression with L1 penalty, (iii) Random Forest Classifier with a
642 maximum depth of 3 and 250 decision stumps (iv) Random Forest with a maximum depth of 6 and 250
643 decision stumps⁸². To treat the problem of class imbalance, a class weight was parameterized according to
644 the percentage of positive (lytic) interactions performed by this phage. The weight of the negative class was
645 always set to 1. For phages with more than 60% of positive interactions, the weight of the positive class
646 was set to 0.8 ; for phages with between 40% and 60% of positive interactions it was set to 1 ; between 30%
647 and 40% of positive interactions positive class weight was set to 1.5 ; between 20% and 30% of positive
648 interactions, class weight was set to 2 and below 20% of positive interactions, class weight was set to 3. All
649 models were provided with the same input features: bacterial core phylogeny and bacterial adsorption
650 factors. Thus a distinct model was used for each individual phage based only on the bacterial adsorption
651 factors and the bacterial core phylogeny (no data regarding defense systems was provided as inputs to the
652 models). Bacterial core phylogeny was embedded using the UMAP algorithm into 8 coordinates⁴³ and
653 adsorption factors were one hot encoded. To ease model training, adsorption factor sharing of each
654 bacterium with the host strain of the phage was also explicitly encoded into features. To avoid model
655 overparameterization, categorical features with less than 3 observations were lumped. This led to a total
656 number of 119 input features to each model. All models were fitted with the scikit-learn python module
657 version 1.1.2. Each fitted model was evaluated using Group 10-Fold Cross-Validation in which train and
658 validation sets were split following bacterial blocks defined by core genome similarity. In this setup,
659 bacterial strains are clustered according to their phylogenetic relationships such that two closely related
660 strains will always be grouped together, either in the train or in the validation set. In our specific case,
bacteria with a core phylogenetic distance of less than 10⁻⁴ substitution per position were grouped together in the same block. This allows us to evaluate the capacity of the model to generalize beyond a certain phylogenetic distance, potentially to completely new ST. Importantly, all the interactions of a given bacterial strain are also always grouped together in the train or in the validation set, forcing the model to learn to predict phage-bacteria interactions based solely on genomic traits (in this specific case, the adsorption factors and the core phylogeny) and not based on previously observed interactions⁸³. The specific value of 10⁻⁴ substitutions/core genome position was chosen as a threshold to cluster together bacterial strains into cross-validation groups for two reasons. First, we previously observed that removing pairs of bacteria below this threshold allowed us to remove the significance of correlation between host range similarity and core phylogenetic distance. We estimated the minimum phylogenetic distance that can be obtained when sampling at random two bacterial strains from different collections (i.e. from different locations, different environments, at different time points by different experimentalists) and observed that 10⁻⁴ substitutions per position was a frequently observed distance and was thus a realistic to measure the capacity of the models to generalize to previously unseen bacterial strains. For each phage, the average AUROC obtained on the validation sets of each Cross-Validation fold were compared and the best model was kept. Model performance was assessed using Precision (True Positives / (True Positives + False Positives)), Recall (True Positives / (True Positives + False Negatives)), F1-score and Average Precision (area under the Precision-Recall curve).

## Building a recommender system for phage cocktails

The recommender system was designed to contain four successive prediction steps. At each step, all phages are ranked based on their predicted probability of infecting the query regarding a step-specific criterion/decision rule and phages with a predicted probability of infection greater than 0.5 are recommended.

The four steps are organized sequentially in order to give priority to the upstream steps compared to the downstream ones as the former were empirically observed to be more precise than the latter. The first step of the pipeline ("Similar as isolation strain") aims at comparing the query strain with the isolations strains in the *Picard* collection based on their adsorption factors and haplotype (ST/O-antigen serotype/H-antigen serotype). Concretely, one machine learning model (Random Forest with a max depth of 4, 80 estimators and a class weight of 1.5 for the positive class) is trained and takes as input features the identity of the phage, as well as boolean features encoding whether the query strain has the same ST/O-antigen/H-antigen/outer core LPS/ABC-dependent capsule serotype as the isolation strain of the phage. The specific decision rules to predict if the query strain is infected by the phage, based on its resemblance with the isolation strain of the phage, were automatically learned by the model at training time on the original interaction matrix. The second step of the pipeline ("Similar as strains in the *Picard* collection") of the pipeline compares the query strain to the rest of the strains in the *Picard* collection based on their haplotypes (combination of ST/O-antigen/H-antigen), using a *k*-neighbors-like approach. For a given query strain, we seek to identify strains in the *Picard* collection which have the same haplotype (same ST, same O-antigen and same H-antigen). If such strains are identified, phages which infect the majority of the bacteria in this group are recommended. The third step ("Trait matching") consists in providing to a random forest classifier as in the first step (max depth of 4; 80 estimators and a class weight of 1.5 for the positive class) all features regarding the bacterial core genome and adsorption factors as well as the identity of the phage. The fourth step of the pipeline is reached if the algorithm does not manage to recommend three phages upstream in the pipeline. A naive (uninformed) strategy is then adopted: the pipeline fills the remaining
704 places in the cocktail with the most generalist phages in the *Guelin* collection (highest number of lytic
705 interactions in the interaction matrix with the Picard collection). The rationale behind this last step is that
706 when the implemented decision rules implemented in the three first steps of the pipeline do not capture the
707 query strain, the only rational decision is to rely on the best uninformed decision rule we have, *i.e.* to
708 recommend phages that lyse the highest number of bacteria in the training set similar to what was done to
709 choose the baseline cocktail. During the process, we explicitly forbid two phages from the same genus and
710 isolated on the same host to be recommended to a given query strain. This design rule guarantees some
711 diversity in the phage cocktail and, in the case of erroneous predictions, avoids recommending three similar
712 that will fail together.

713 **Statistics and reproducibility**

714 Statistical tests were performed using either python packages scipy.stats v1.11.2<sup>84</sup>, statsmodels v0.14.1 or
715 in R using VEGAN<sup>38</sup> v2.6-4, MCMCglmm<sup>45</sup> v2.35.

716 No statistical method was used to predetermine sample size. No data were excluded from the analyses. The
717 experiments were not randomized. The investigators were not blinded when assessing phage-bacteria
718 interactions but were unaware of the identity of phages and bacteria when annotating the dataset. All data
719 points were treated and analyzed equally.

720 **Citation diversity**

721 To promote gender equality and inclusivity in research, we are convinced of the importance of
722 acknowledging gender bias in research article citation. Using a custom script available on our github
723 (https://github.com/mdmparis/Estimating gender bias in references), we estimated that among the 83
724 references cited in the main text, approximately 30% (25) have a female first author and approximately
725 14% (12) have a female last author.

726 **Data availability**

727 All the data used to generate the analyses presented here is available on the Github repository of the project
728 (https://github.com/mdmparis/coli_phage_interactions_2023). Raw images of the plaque assay
729 experiments are available on the Zenodo repository of the project (https://zenodo.org/records/10202713).
730 Bacterial genomes are available on the Figshare repository under the DOI 10.6084/m9 figshare.25941691.
731 The encoded interaction matrix is also available in the Viral Host Range database
732 (https://viralhostrangedb.pasteur.cloud/).

733 **Code availability**

734 All the code used to generate the analyses is available on the GitHub repository of the projects
735 (https://github.com/mdmparis/coli_phage_interactions_2023).

736 **Acknowledgements**

737 We are grateful to members of the MDM lab, Beatriz Beamud, David Bikard, Ariel Lindner and Vincent
738 Libis for their useful comments on earlier versions of the manuscript. We thank Damien Maura, Mathieu
739 Galtier, Marine Henry, Nicolas Dufour, Raphaëlle Delattre and Marta Lourenço from the Bacteriophage,
740 Bacterium, Host lab of Institut Pasteur for the isolation of phages of the *Guelin* collection. We thank the
741 ColoColi group for sharing the strains of the *ColocColi* collection. We thank Juliette Meyer from the
742 Bioinformatic hub of Institut Pasteur for advice on statistical analysis. Several bioinformatic analyses were
743 performed on the Core Cluster of the Institut Français de Bioinformatique (IFB) (ANR-11-INBS-0013).
744 We are also grateful to Sandra Legout and the team managing scientific archives at the Institut Pasteur for
745 their kind help in naming the phage collection.

746 Author Contributions Statement

747 B.G. and A.B. conceptualized the study and secured funding. B.G., J.B., T.D., O.C. collected, checked,
748 and stored phages and strains collections. B.G., I.C., N.D., J.B. performed phage-bacteria interactions
749 experiments. B.G., I.C., F.T., H.G., O.C. performed sequencing experiments. F.T. performed the
750 bioinformatics analysis of the genomes. H.V. performed statistical inferences and trained the machine
751 learning models. B.G., H.V., F.T. analyzed and interpreted the data. A.B., B.G., H.V., F.T. wrote the first
752 draft. All authors reviewed and edited the final draft.

753 B. G., H. V., F. T., I. C., N. D. and A.B. are supported by the CRI Research Fellowship to A.B. from the
754 Bettencourt Schueller Foundation, the ATIP-Avenir program from INSERM (R21042KS/RSE22002KSA),
755 the Emergence program from the University of Paris-Cité (RSFVJ21IDXB6_DANA) ERC Starting Grant
756 (PECAN 101040529) and the core funding of Institut Pasteur. E.D. was partially supported by the
757 “Fondation pour la Recherche Médicale” (Equipe FRM 2016, grant number DEQ20161136698). L.D. was
758 partially supported by Agence Nationale de la Recherche (ANR-19-AMRB-0002 and ANR-20-CE92-
759 0048). B.G. was partially supported by Agence Nationale de la Recherche (ANR-19-AMRB-0002).

760 Competing Interests Statement

761 The authors declare no competing interests.
762 **References**

763 1. Kortright, K. E., Chan, B. K., Koff, J. L. & Turner, P. E. Phage Therapy: A Renewed Approach to Combat Antibiotic-Resistant Bacteria. *Cell Host Microbe* **25**, 219–232 (2019).

765 2. Strathdee, S. A., Hatfull, G. F., Mutalik, V. K. & Schooley, R. T. Phage therapy: From biological mechanisms to future directions. *Cell* **186**, 17–31 (2023).

767 3. Lood, C. et al. Digital phagograms: predicting phage infectivity through a multilayer machine learning approach. *Curr. Opin. Virol.* **52**, 174–181 (2022).

769 4. Nobrega, F. L. et al. Targeting mechanisms of tailed bacteriophages. *Nat. Rev. Microbiol.* **16**, 760–773 (2018).

771 5. Georjon, H. & Bernheim, A. The highly diverse antiphage defence systems of bacteria. *Nat. Rev. Microbiol.* **21**, 686–700 (2023).

773 6. Maffei, E. et al. Systematic exploration of Escherichia coli phage-host interactions with the BASEL phage collection. *PLoS Biol.* **19**, e3001424 (2021).

775 7. Stanley, S. Y. & Maxwell, K. L. Phage--encoded Anti-CRISPR Defenses. *Annu. Rev. Genet.* **52**, 445–464 (2018).

777 8. Krüger, D. H. & Bickle, T. A. Bacteriophage survival: multiple mechanisms for avoiding the deoxyribonucleic acid restriction systems of their hosts. *Microbiol. Rev.* **47**, 345–360 (1983).

779 9. Bertozzi Silva, J., Storms, Z. & Sauvageau, D. Host receptors for bacteriophage adsorption. *FEMS Microbiol. Lett.* **363**, fnw002 (2016).

781 10. Tesson, F. et al. Systematic and quantitative view of the antiviral arsenal of prokaryotes. *Nat. Commun.* **13**, 2561 (2022).

783 11. Piel, D. et al. Phage-host coevolution in natural populations. *Nat. Microbiol.* **7**, 1075–1086 (2022).

784 12. Kauffman, K. M. et al. Resolving the structure of phage-bacteria interactions in the context of natural diversity. *Nat. Commun.* **13**, 372 (2022).

786 13. Korf, I. H. E. et al. Still Something to Discover: Novel Insights into Escherichia coli Phage Diversity and Taxonomy. *Viruses* **11**, 454 (2019).

788 14. Walsh, S. K. et al. The host phylogeny determines viral infectivity and replication across Staphylococcus host species. *PLoS Pathog.* **19**, e1011433 (2023).

790 15. Cuervo, A. et al. Structures of T7 bacteriophage portal and tail suggest a viral DNA retention and ejection mechanism. *Nat. Commun.* **10**, 3746 (2019).

792 16. Hu, B., Margolin, W., Molineux, I. J. & Liu, J. The bacteriophage t7 virion undergoes extensive structural remodeling during infection. *Science* **339**, 576–579 (2013).

794 17. Suga, A., Kawaguchi, M., Yonesaki, T. & Otsuka, Y. Manipulating Interactions between T4 Phage Long Tail Fibers and Escherichia coli Receptors. *Appl. Environ. Microbiol.* **87**, e0042321 (2021).

796 18. Srikant, S., Guegler, C. K. & Laub, M. T. The evolution of a counter-defense mechanism in a virus constrains its host range. *eLife* **11**, e79549 (2022).

798 19. Abby, S. S., Néron, B., Ménager, H., Touchon, M. & Rocha, E. P. C. MacSyFinder: a program to mine genomes for molecular systems with an application to CRISPR-Cas systems. *PLoS One* **9**, e110726 (2014).

801 20. Néron, B. et al. MacSyFinder v2: Improved modelling and search engine to identify molecular systems in genomes. *Peer Community J.* **3**, (2023).

803 21. Boeckaerts, D., Stock, M., De Baets, B. & Briers, Y. Identification of Phage Receptor-Binding Protein Sequences with Hidden Markov Models and an Extreme Gradient Boosting Classifier. *Viruses* **14**, 1329 (2022).

806 22. Pan, J. et al. GSPHI: A novel deep learning model for predicting phage-host interactions via multiple biological information. *Comput. Struct. Biotechnol. J.* **21**, 3404–3413 (2023).

808 23. Wang, Y. et al. An Effective Model for Predicting Phage-Host Interactions Via Graph Embedding Representation Learning With Multi-Head Attention Mechanism. *IEEE J. Biomed. Health Inform.* **27**, 3061–3071 (2023).

811 24. Boeckaerts, D. et al. Predicting bacteriophage hosts based on sequences of annotated receptor-binding proteins. *Sci. Rep.* **11**, 1467 (2021).
813 25. Bajiya, N., Dhall, A., Aggarwal, S. & Raghava, G. P. S. Advances in the field of phage-based therapy with special emphasis on computational resources. *Brief. Bioinform.* **24**, bbac574 (2023).

815 26. Moller, A. G. et al. Genes Influencing Phage Host Range in Staphylococcus aureus on a Species-Wide Scale. *mSphere* **6**, e01263-20 (2021).

817 27. Beamud, B. et al. Genetic determinants of host tropism in Klebsiella phages. *Cell Rep.* **42**, 112048 (2023).

819 28. Haudiquet, M., Buffet, A., Rendueles, O. & Rocha, E. P. C. Interplay between the cell envelope and mobile genetic elements shapes gene flow in populations of the nosocomial pathogen Klebsiella pneumoniae. *PLoS Biol.* **19**, e3001276 (2021).

822 29. Boeckaerts, D. et al. Prediction of Klebsiella phage-host specificity at the strain level. *Nat. Commun.* **15**, 4355 (2024).

824 30. Keith, M. et al. Predictive phage therapy for Escherichia coli urinary tract infections: cocktail selection for therapy based on machine learning models. 2023.11.23.568453 Preprint at https://doi.org/10.1101/2023.11.23.568453 (2023).

827 31. Antimicrobial Resistance Collaborators. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. *Lancet Lond. Engl.* **399**, 629–655 (2022).

829 32. Tenaillon, O., Skurnik, D., Picard, B. & Denamur, E. The population genetics of commensal Escherichia coli. *Nat. Rev. Microbiol.* **8**, 207–217 (2010).

831 33. Denamur, E., Clermont, O., Bonacorsi, S. & Gordon, D. The population genetics of pathogenic Escherichia coli. *Nat. Rev. Microbiol.* **19**, 37–54 (2021).

833 34. Galardini, M. et al. Phenotype inference in an Escherichia coli strain panel. *eLife* **6**, e31035 (2017).

834 35. Galardini, M. et al. Major role of iron uptake systems in the intrinsic extra-intestinal virulence of the genus Escherichia revealed by a genome-wide association study. *PLoS Genet.* **16**, e1009065 (2020).

836 36. Bolduc, B. et al. vConTACT: an iVirus tool to classify double-stranded DNA viruses that infect Archaea and Bacteria. *PeerJ* **5**, e3243 (2017).

838 37. O'Leary, N. A. et al. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. *Nucleic Acids Res.* **44**, D733-745 (2016).

840 38. Dixon, P. VEGAN, a package of R functions for community ecology. *J. Veg. Sci.* **14**, 927–930 (2003).

842 39. Smug, B. J., Szczepaniak, K., Rocha, E. P. C., Dunin-Horkawicz, S. & Mostowy, R. J. Ongoing shuffling of protein fragments diversifies core viral functions linked to interactions with bacterial hosts. 2022.12.27.521992 Preprint at https://doi.org/10.1101/2022.12.27.521992 (2023).

845 40. Pas, C., Latka, A., Fieseler, L. & Briers, Y. Phage tailspike modularity and horizontal gene transfer reveals specificity towards E. coli O-antigen serogroups. *Virol. J.* **20**, 174 (2023).

847 41. Sørensen, A. N., Woudstra, C., Sørensen, M. C. H. & Brøndsted, L. Subtypes of tail spike proteins predicts the host range of Ackermannviridae phages. *Comput. Struct. Biotechnol. J.* **19**, 4854–4867 (2021).

850 42. Gencay, Y. E., Gambino, M., Prüssing, T. F. & Brøndsted, L. The genera of bacteriophages and their receptors are the major determinants of host range. *Environ. Microbiol.* **21**, 2095–2111 (2019).

852 43. McInnes, L., Healy, J., Saul, N. & Großberger, L. UMAP: Uniform Manifold Approximation and Projection. *J. Open Source Softw.* **3**, 861 (2018).

854 44. Gilchrist, C. L. M. & Chooi, Y.-H. clinker & clustermap.js: automatic generation of gene cluster comparison figures. *Bioinforma. Oxf. Engl.* **37**, 2473–2475 (2021).

856 45. Hadfield, J. D. MCMC Methods for Multi-Response Generalized Linear Mixed Models: The MCMCglmm R Package. *J. Stat. Softw.* **33**, (2010).

858 46. Casjens, S. R. & Hendrix, R. W. Bacteriophage lambda: early pioneer and still relevant. *Virology* **0**, 310–330 (2015).

860 47. Guo, X., Yin, Y., Dong, C., Yang, G. & Zhou, G. On the Class Imbalance Problem. in *2008 Fourth International Conference on Natural Computation* vol. 4 192–201 (2008).
862 48. Abedon, S. T., Danis-Wlodarczyk, K. M. & Wozniak, D. J. Phage Cocktail Development for Bacteriophage Therapy: Toward Improving Spectrum of Activity Breadth and Depth. *Pharm. Basel* **14**, 1019 (2021).

863
864
865 49. Holtappels, D., Alfenas-Zerbini, P. & Koskella, B. Drivers and consequences of bacteriophage host range. *FEMS Microbiol. Rev.* **47**, fuad038 (2023).

867 50. Rocha, E. P. C. & Bikard, D. Microbial defenses against mobile genetic elements and viruses: Who defends whom from what? *PLOS Biol.* **20**, e3001514 (2022).

869 51. Rousset, F. et al. Phages and their satellites encode hotspots of antiviral systems. *Cell Host Microbe* **30**, 740-753.e5 (2022).

871 52. La Combe, B. et al. Pneumonia-Specific Escherichia coli with Distinct Phylogenetic and Virulence Profiles, France, 2012-2014. *Emerg. Infect. Dis.* **25**, 710–718 (2019).

873 53. Debarbieux, L. et al. Bacteriophages can treat and prevent Pseudomonas aeruginosa lung infections. *J. Infect. Dis.* **201**, 1096–1104 (2010).

875 54. Maura, D. et al. Intestinal colonization by enteroaggregative Escherichia coli supports long-term bacteriophage replication in mice. *Environ. Microbiol.* **14**, 1844–1854 (2012).

877 55. Dufour, N., Debarbieux, L., Fromentin, M. & Ricard, J.-D. Treatment of Highly Virulent Extraintestinal Pathogenic Escherichia coli Pneumonia With Bacteriophages. *Crit. Care Med.* **43**, e190-198 (2015).

880 56. Dufour, N. et al. Bacteriophage LM33_P1, a fast-acting weapon against the pandemic ST131-O25b:H4 Escherichia coli clonal complex. *J. Antimicrob. Chemother.* **71**, 3072–3080 (2016).

882 57. Galtier, M. et al. Bacteriophages to reduce gut carriage of antibiotic resistant uropathogens with low impact on microbiota composition. *Environ. Microbiol.* **18**, 2237–2245 (2016).

884 58. Galtier, M. et al. Bacteriophages Targeting Adherent Invasive Escherichia coli Strains as a Promising New Treatment for Crohn's Disease. *J. Crohns Colitis* **11**, 840–847 (2017).

886 59. De Sordi, L., Khanna, V. & Debarbieux, L. The Gut Microbiota Facilitates Drifts in the Genetic Diversity and Infectivity of Bacterial Viruses. *Cell Host Microbe* **22**, 801-808.e3 (2017).

888 60. Lourenço, M. et al. The Spatial Heterogeneity of the Gut Limits Predation and Fosters Coexistence of Bacteria and Bacteriophages. *Cell Host Microbe* **28**, 390-401.e5 (2020).

890 61. Chiu, C.-L., Clack, N. & Community, T. N. napari: a Python Multi-Dimensional Image Viewer Platform for the Research Community. *Microsc. Microanal.* **28**, 1576–1577 (2022).

892 62. Lamy-Besnier, Q., Brancotte, B., Ménager, H. & Debarbieux, L. Viral Host Range database, an online tool for recording, analyzing and disseminating virus-host interactions. *Bioinformatics* **37**, 2798–2801 (2021).

895 63. Trivedi, U. H. et al. Quality control of next-generation sequencing data without a reference. *Front. Genet.* **5**, 111 (2014).

897 64. Bankevich, A. et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. *J. Comput. Biol. J. Comput. Mol. Cell Biol.* **19**, 455–477 (2012).

899 65. Seemann, T. Prokka: rapid prokaryotic genome annotation. *Bioinforma. Oxf. Engl.* **30**, 2068–2069 (2014).

901 66. Jolley, K. A., Bray, J. E. & Maiden, M. C. J. Open-access bacterial population genomics: BIGSdb software, the PubMLST.org website and their applications. *Wellcome Open Res.* **3**, 124 (2018).

903 67. Beghain, J., Bridier-Nahmias, A., Le Nagard, H., Denamur, E. & Clermont, O. ClermonTyping: an easy-to-use and accurate in silico method for Escherichia genus strain phylotyping. *Microb. Genomics* **4**, e000192 (2018).

906 68. Perrin, E. P. C., A., . Rocha. PanACoTA: a modular tool for massive microbial comparative genomics. *NAR Genomics Bioinforma* **3**, (2021).

908 69. Letunic, I. & Bork, P. Interactive Tree Of Life (iTOL) v5: an online tool for phylogenetic tree display and annotation. *Nucleic Acids Res.* **49**, W293–W296 (2021).

910 70. Moraru, C., Varsani, A. & Kropinski, A. M. VIRIDIC-A Novel Tool to Calculate the Intergenomic Similarities of Prokaryote-Infecting Viruses. *Viruses* **12**, 1268 (2020).
912 71. Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular
913 interaction networks. Genome Res. **13**, 2498–2504 (2003).

914 72. Bessonov, K. et al. ECTyper: in silico Escherichia coli serotype and species prediction from raw and
915 assembled whole-genome sequence data. *Microb. Genomics* **7**, 000728 (2021).

916 73. Rendueles, O., Garcia-Garcerà, M., Néron, B., Touchon, M. & Rocha, E. P. C. Abundance and co-
917 occurrence of extracellular capsules increase environmental breadth: Implications for the emergence
918 of pathogens. *PLoS Pathog.* **13**, e1006525 (2017).

919 74. Mirdita, M., Steinegger, M. & Söding, J. MMseqs2 desktop and local web server app for fast,
920 interactive sequence searches. *Bioinforma. Oxf. Engl.* **35**, 2856–2858 (2019).

921 75. Lam, M. M. C., Wick, R. R., Judd, L. M., Holt, K. E. & Wyres, K. L. Kaptive 2.0: updated capsule
922 and lipopolysaccharide locus typing for the Klebsiella pneumoniae species complex. *Microb.
923 Genomics* **8**, 000800 (2022).

924 76. Rodríguez-Gironés, M. A. & Santamaría, L. A new algorithm to calculate the nestedness temperature
925 of presence–absence matrices. *J. Biogeogr.* **33**, 924–935 (2006).

926 77. Barber, M. J. Modularity and community detection in bipartite networks. *Phys. Rev. E* **76**, 066102
927 (2007).

928 78. Eddy, S. R. Accelerated Profile HMM Searches. *PLOS Comput. Biol.* **7**, e1002195 (2011).

929 79. Mistry, J. et al. Pfam: The protein families database in 2021. *Nucleic Acids Res.* **49**, D412–D419
930 (2021).

931 80. Lin, Z. et al. Evolutionary-scale prediction of atomic-level protein structure with a language model.
932 *Science* **379**, 1123–1130 (2023).

933 81. van Kempen, M. et al. Fast and accurate protein structure search with Foldseek. *Nat. Biotechnol.* **1–4**
934 (2023) doi:10.1038/s41587-023-01773-0.

935 82. Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. *Mach. Learn. PYTHON*.

936 83. Stock, M., Pahikkala, T., Airola, A., Waegeman, W. & De Baets, B. Algebraic shortcuts for leave-
937 one-out cross-validation in supervised network inference. *Brief. Bioinform.* **21**, 262–271 (2020)

938 84. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat
939 Methods **17**, 261–272 (2020)
# Supplementary Materials

## Table of contents:

### Supplementary Tables

- Supplementary Table 1: *Picard* Collection

- Supplementary Table 2: *Guelin* Collection

- Supplementary Table 3: Interaction matrix

- Supplementary Table 4: Outer membrane proteins of strains from the *Picard* collection

- Supplementary Table 5: Defense systems of strains from the *Picard* collection

- Supplementary Table 6: Test collection

- Supplementary Table 7: Phage cocktails recommended for the test collection

### Supplementary Texts

- Supplementary Text 1 – Short biography of Bertrand Picard

- Supplementary Text 2 - Bertrand Picard collection encompasses the *Escherichia* genus diversity

- Supplementary Text 3 – Short biography of Antonina Guelin

- Supplementary Text 4 – Link between phage genome size and number of lytic interactions within phage morphotypes

- Supplementary Text 5 – Nestedness and modularity properties of the interaction matrix

- Supplementary Text 6 – RBP sharing is sufficient to explain host range similarity among phages, even among different genera

- Supplementary Text 7 – Link between the number of lytic interactions performed by a phage and the number of RBPs encoded in its genome

- Supplementary Text 8 – Multivariate analysis of the bacterial traits associated with lytic interactions in the matrix

- Supplementary Text 9 – Link between bacterial susceptibility to phages and bacterial genome size or number of antimicrobial genes encoded in the genome

- Supplementary Text 10 – Rationale for the number of phages used in our phage cocktail recommender system

- Supplementary Text 11 – Assessment of antagonism or synergy between phages in cocktails

- Supplementary Text 12 - Limit of genomic detection

### Supplementary Figures

- Supplementary Figure 1. Interactions between the three laboratory-adapted strains MG1655, BE, and BL21 with the phages of the *Antonina* Guelin collection.
# Supplementary Texts

## Supplementary Text 1 – Short biography of Bertrand Picard

Born in 1950, Bertrand Picard is a French medical microbiologist who devoted most of his career to studying *Escherichia coli* natural isolate diversity and the link between strain genetic features and lifestyle. He obtained his MD degree from the University Paris VII in 1978, followed by a PhD in Biochemistry and Molecular Biology in 1990 at the University Paris V. He completed the molecular taxonomy course at the Institut Pasteur in 1988. He spent most of his hospital career as a medical bacteriologist at the Assistance Publique-Hôpitaux de Paris as an Assistant (1980), then Associate (1986), and finally Full Professor (1994), with a 10-year period in the University hospitals of Brest in Britany.

His research was based on two main intricate objectives. First, he characterized the structural basis of the polymorphism of the carboxylesterase B of *E. coli*. Second, by studying numerous natural isolates from different contexts (commensal, pathogen), in different hosts (humans, non-human mammals, birds) from different countries, he looked for genetic variants associated with a specific lifestyle. His major contribution was to show that strains exhibiting a carboxylesterase B₂ electrophoretic type were more frequently found in extraintestinal infections, exhibited numerous virulence-associated genes, and were highly virulent in a mouse model of sepsis. He is the author of over 150 publications in peer-reviewed journals.

## Supplementary Text 2 - Bertrand Picard collection encompasses the *Escherichia* genus diversity

*Escherichia phylogenetic diversity.* We gathered strains from six species and eight phylogroups of the *Escherichia* genus (*Escherichia coli*, *n* = 378). To estimate genomic diversity in a finer granularity, we also analyzed Sequence Types (ST). In the pre-genomic era, the clonal population structure of the *Escherichia* genus was typically described in terms of ST, which was determined by MultiLocus Sequence Typing (MLST)¹,². *In silico* MLST revealed that our collection accounts for 163 STs, with 50% of the collection covered by 14 main STs (see **Figure 1A** and **1B**)³. Beyond ST, we analyzed additional loci known to frequently recombine, such as the *rfb* locus responsible for the production of the O-antigen and the *fli* locus responsible for the production of the H-antigen. Both O- and H-antigens are surface polysaccharides that undergo strong diversifying selection and are involved in phage susceptibility¹,⁴. O and H-antigen serotypes were predicted *in silico*³. A total of 93 O-antigen (out of approximately 180 currently described in *E. coli*) and 41 H-antigen serotypes (out of approximately 50) are present in the *Picard* collection, including all the most abundant ones (see **Figure 1C**)⁵,⁶. Apart from phylogenetic diversity, we paid attention to the pathotype (commensal, environmental, pathogenic) (see **Figure 1D**) and ecological source diversities of the isolates (see **Figure 1E**).

*Escherichia genomic diversity of phage interaction traits.* Using the genomic sequence of the strains in the *Picard* collection, we sought to characterize the traits that could potentially influence phage susceptibility. We identified the bacterial surface structures known to be potential receptors for phages⁷,⁸. In *Escherichia*, many surface polysaccharides were shown to be involved in host cell recognition by phages. Out of these, the lipopolysaccharide (LPS) is a highly diverse glycolipid covering more than 66% of the bacterial outer surface. It is composed of the lipid-A, the LPS inner and outer cores, and the highly diverse O-antigen described earlier⁹. The latter two are frequently described as phage receptors in *E. coli*⁸. In our collection, we identified the five major groups of LPS outer core (see **Figure 1F**). Capsules have also been previously described as receptors for capsule serotype-dependent phages in *E. coli*¹⁰. We developed an approach to identify strains encoding an ABC-dependent capsule and to cluster them according to their K-serotype (see methods). We observed that most of the main K-serotypes described as phage receptors were well
represented in the *Picard* collection (see **Figure 1G**). Capsules are well known to shape the phage host range in the *K. pneumoniae* species<sup>11,12</sup> and *Klebsiella* capsules have been previously identified in *E. coli* natural isolates<sup>13</sup>. We hypothesized that these capsules could also influence phage susceptibility in *E. coli* and used Kaptive to detect the *Klebsiella* capsule-encoding strains in our collection<sup>14</sup>. We found that 22 of our 403 *Escherichia* isolates (5.4%) encode such capsules (strain serotypes available in **Supplementary Table 1**). All the *K. pneumoniae* capsule encoding strains belonged to phylogroups A, B1, or C and had O-antigen serotypes O8, O9, or O89, in perfect agreement with a preceding study on the distribution of *Klebsiella* capsules among environmental *E. coli* isolates (see **Supplementary Table 1**)<sup>13</sup>. We also detected twelve genes coding for the outer membrane proteins that were recurrently identified in the literature as putative phage receptors for *E. coli* (see **Supplementary Table 4**)<sup>8</sup>, revealing that they are ubiquitous in *Escherichia* genomes (each one was found in more than 97% of the collection). Phage adsorption can, however, be easily averted by a single substitution in the gene coding for the phage receptor. When looking for the genomic diversity of genes coding for these proteins, we thus identified a large diversity of genomic variants (see **Figure 1H**). Gene variants were clustered at 99% genomic identity, meaning that two variants within the same cluster diverged by at most approximately 5 amino acids from each other. Altogether, these analyses provide an extensive catalog of combinations of putative phage adsorption factors within the *Picard* collection. We next focused on intra-cellular defense. Using DefenseFinder<sup>15</sup>, we identified known antiphage systems (see **Supplementary Table 5**). On average, an isolate encodes 8 defense system subfamilies, with considerable heterogeneity among strains (i.e. from 1 to 16 defense systems) (see **Figure 1I**). The pan immune system of the *Escherichia* genus comprises a total of 137 antiphage defense systems. Some of them are widespread (e.g. MazEF or RM type I), while many are present in only a few isolates (see **Figure 1J**). As a result, our genomic analysis of the *Picard* collection of *Escherichia* strains provides a large panel of genomic traits relevant to investigating phage-bacteria interactions at the scale of the genus.

**Supplementary Text 3 – Short biography of Antonina Guelin**

Antonina Guelin (born Chtdedrina) was a French medical doctor and biologist of Russian origin. Born in 1904, in Saint-Petersburg (Russia), she died in 1988. She was first Assistant at the Embryology Laboratory of the Collège de France (1931-1934), where her work focused on the chicken embryo. She then entered the laboratory in the anaerobes department (1934-1941), and the bacteriophage department (1941-1968) headed by Eugène Wollman at the Pasteur Institute. Due to her early works on anaerobic strains, she pioneered the research on phage of *Clostridium perfringens* and *C. histolyticum*. During this period, she notably demonstrated the benefits of monitoring bacteriophages to assess the sanitary quality of fresh water, and the role of bacteriophages in the regulation of bacterial population in aquatic environments. For those works, she received an award from the French National Academy of Medicine. She finished her career as a researcher at the Roscoff Biological Station in Britany (France) (1970-1988).

**Supplementary Text 4 – Link between phage genome size and number of lytic interactions within phage morphotypes**

We wondered if the tendency of Myovirus to have a broader host range may be due to their larger genome allowing them to adapt to more diverse host genetic backgrounds at any step of their life cycle (adsorption to host, intracellular replication, etc.). However, when taking the phage morphotype into account, phage genome size did not show a significant association with the number of lytic interactions performed by each phage (Linear model: genome size = 2.3 x 10⁻⁴, *p*=0.17).
**Supplementary Text 5 – Nestedness and modularity properties of the interaction matrix**

A frequently discussed macroscopic property of bipartite networks like phage-bacteria interaction matrices is how the interactions are distributed within the network: from modular (independent groups of phages infect independent groups of bacteria) to nested (host range of more specific phages/bacteria tend to be subsets of the host range of more generalist ones)¹⁶. Our interaction matrix exhibited strong nestedness (BINMATNEST algorithm, T=11.63, *p*<10⁻⁵) and weak modularity (BRIM algorithm, Q=0.23, *p*<10⁻³) (Supplementary Figure 2).

**Supplementary Text 6 – RBP sharing is sufficient to explain host range similarity among phages, even among different genera**

We hypothesized that RBP could be the major phage trait explaining the infection patterns observed in the interaction matrix. To explore this hypothesis, we identified a triplet of phages that exhibit very similar infection patterns while belonging to different genera and being isolated on three distinct isolation strains (Supplementary Figure 4A). We sought to identify the genetic factors underlying the similarity in their infection patterns and wondered if this similarity could be explained by shared RBP. Two of these phages (T145_P2 and AL505_P3) belong to the *Vectrevirus* genus (*Autographiviridae/Molineuxvirinae*), but the third one (MT1B1_P3) is a *Kayfunavirus* (*Autographiviridae/Studiervirinae*) and shared very few genomic similarities with the other ones. The three phages shared only three of their core genes: (i) their large terminase subunit, (ii) their RNA polymerase, which is also conserved among all *Vectreviruses* and *Kayfunaviruses*, and (iii) a gene identified as an RBP in the *Kayfunavirus* MT1B1_P3 and in both *Vectreviruses*. While sharing limited amino acid identity (56% between the *Vectrevirus* and *Kayfunavirus* cognate gene), domain annotation using the PFAM-A database and the structure prediction of each RBP monomer allowed us to identify that the three RBP of interest shares a common architecture (Supplementary Figure 4B)¹⁷,¹⁸. Structural homology search revealed that this protein architecture can be found in other phage tailspikes, such as in the K1F phage tailspike (K1F TSP), which has been shown to use K1-serotype capsule as a target¹⁹⁻²¹ (Supplementary Figure 4B). K1F TSP has been shown to bind to polysialic acid and have endosialidase activity²¹. We hypothesized that such RBP could, therefore, also use K1 capsules as a receptor. We identified 45 strains encoding a K1-serotype capsule in the *Picard* collection, including 4 isolation strains, 3 of which were precisely the isolation strains of phages T145_P2, AL505_P3, and MT1B1_P3. When testing the decision rule that our three phages of interest infect exactly the K1-encapsulated strains, we obtained a precision (percentage of true positives among all predicted lytic interactions) of 73% and a recall (percentage of true positives among all lytic interactions) of 72%. As a comparison, we applied the same decision rule – lytic interaction of a bacterial strain if and only if the strain harbors a K1-serotype capsule – to all the other *Vectreviruses* and *Kayfunaviruses* and obtained a precision of 7% and a recall of 14%. This suggests that T145_P2, MT1B1_P3, and AL55_P3 exhibit important similarities in their infection patterns because they possess a similar RBP in spite of being isolated and propagated on distinct bacterial strains. This detailed example further shows that phages that share their RBP, regardless of their isolation strain, also have similar infection patterns. Altogether, this further strengthens the idea that RBPs are the major determinant of phage host range.

**Supplementary Text 7 – Link between the number of lytic interactions performed by a phage and the number of RBPs encoded in its genome**
The systematic identification of RBP in all our phages also allowed us to wonder if the number of lytic
interactions performed by a phage could be linked with the number of RBPs it encodes. We thus tested
whether the number of lytic interactions performed by a phage host range would be the consequence of the
number of RBP it encodes in its genome. Under this hypothesis, we would expect that phages that perform
more lytic interactions (e.g. Myoviruses, which perform up to 40% of lytic interactions on the *Picard*
collection) would have more RBPs than phages that are specific (e.g. Podo- and Siphovirus which
sometimes perform less than 10% of lytic interactions). We ran a systematic analysis on the whole *Guelin*
collection of phages and observed that all phages in the collection encode either one or two predicted RBP
(See methods)²². This estimate was supported by the analysis of the literature in which we found no
examples of *Caudoviricetes* with more than two RBP-encoding genes, with the notable exception of phages
from the *Ackermannviridae* family harboring four tailspikes in star-like structures²³. None of these phages
were included in the *Guelin* collection.

As a consequence, we did not find any significant correlation between the number of RBP of a phage and
the number of lytic interactions it performs (Linear Mixed Model (LMM), β□□□□□□ = 28.11, 95% credible
interval = [-33.09, 83.58], pMCMC = 0.33). This implies that the number of lytic interactions a phage
performs is primarily influenced by the tropism of its RBP rather than the number of RBP it encodes.

**Supplementary Text 8 – Multivariate analysis of the bacterial traits associated with lytic interactions in the matrix**

To assess the statistical association between bacterial traits and lytic interactions in the matrix, we also
followed a multivariate approach instead of treating each phage infection separately. We fitted a single
model on all the phage-bacteria interactions in the matrix. We first fitted a PERMANOVA model to which
we provided the presence/absence of both defense system subfamilies and adsorption factors as covariates
("full" model). We compared this full model to (i) a submodel having only adsorption factors as covariates
and (ii) a submodel having only defense system subfamilies as covariates. Model comparison was
performed using the Akaike Information Criterion (AIC) (See methods). The best possible fit was obtained
with the adsorption factors-based model (PERMANOVA, AIC=-435, adjusted R²=0.27) compared with the
full model (AIC=-418, adjusted R²=0.28) or the defense systems-based model (AIC=-372, adjusted
R²=0.05). This approach also suggests that adsorption factors better explain bacterial host range than
bacterial defense mechanisms.

**Supplementary Text 9 – Link between bacterial susceptibility to phages and bacterial genome size or number of antimicrobial genes encoded in the genome**

We observed a significant but low correlation between a larger defense arsenal and decreased MLC, given
phage infection (Figure 4D). This correlation was robust to the removal of extreme values and to the
random subsetting of bacterial ST, indicating that it was not driven by a limited number of bacterial clades.
To control for potential confounding factors that would be associated with antiphage systems, we also
assessed the statistical association between the average MLC given infection and (i) the number of
antimicrobial resistance (AMR) genes encoded in each bacterial genome or (ii) the bacterial genome size.
AMR genes were chosen as a control variable as they undergo evolutionary dynamics similar to the ones
of defense systems. Genome size was chosen as it is related to the composition of the pangenome in a given
strain²⁴. None of these variables were *a priori* expected to be associated with susceptibility to phages, but
if they were, this would suggest that the observed correlation between the MLC-given phage infection and
the number of defense systems per genome is more likely to be a “pangenome-related” effect than a true
effect of the presence of defense systems (Supplementary Figure 6). Only genome size was weakly but significantly associated with the average MLC given infection (Genome size: $\beta_{genome size} = -4.49 \times 10^{-7}$, 95% credible interval = [-8.5x10⁻⁷, -1.0x10⁻⁷], pMCMC = 0.05; AMR genes: $\beta_{Num AMR genes} = -0.03$, 95% credible interval = [-0.05, 0.01], pMCMC = 0.11). However, this correlation between the average MLC given infection and the bacterial genome size did not hold once the extreme values in terms of genome size were removed from the dataset. Furthermore, when putting both genome size and the number of defense systems as covariates in a model, only the number of defense systems was significantly associated with the MLC-given phage infection. This indicates that the information carried by the genome size variable was recapitulated by the number of defense systems and that the best fit is obtained with this variable.

Both controls strengthened the idea that the negative correlation between the number of defense systems encoded by a strain and its average MLC given infection is not spurious. It is important to note that although significant, the effect size of the number of defense systems on the MLC-given phage infection was low ($\text{Num DS} = -0.04$).

**Supplementary Text 10 – Rationale for the number of phages used in our phage cocktail recommender system**

Before designing a machine learning-based phage cocktail recommender system, we first set out to determine what should be the size, i.e. the number of phages, in any cocktail. This parameter should be large enough to cover diverse bacteria and provide more robustness to erroneous phage recommendations but remain low to avoid undesired side effects and limit the cost of production. We evaluated the number of infected bacteria by the most covering k-uplet of phages. We observed that the most covering triplet of phages infect 63% of the Picard collection, whereas the most covering 4-uplet of phages only reached 68% of infected bacteria (Supplementary Figure 7). This small increase in efficiency led us to choose a cocktail size of k=3 phages, offering a good trade-off between efficiency and parsimony.

**Supplementary Text 11 – Assessment of antagonism or synergy between phages in cocktails**

We have evaluated the tailored, generic, and baseline cocktails experimentally by plaque assay experiments (see methods). We tested each possible decomposition of the cocktails (each phage individually, each pair of phages, and the full 3-phage cocktail) to identify putative synergistic or antagonistic behavior between phages. Interestingly, no antagonistic nor synergistic phage-phage interaction was observed in this experimental setup.

**Supplementary Text 12**

Our collections, while large, only represent a subset of the diversity of all *Escherichia* natural isolates. Some adsorption factors were not fully identified in our research, including the Enterobacterial Common Antigen (ECA) or the Wzy-dependent capsular serotypes in *E. coli*⁸,²⁵. Genomic analyses also have inherent limitations, as they do not consider gene regulation and expression in response to the bacterial environment. We evaluated phage-bacteria interactions under controlled conditions, specifically in a non-limiting limited medium. The expression of many proteins, including certain outer membrane proteins or capsular components, may exhibit variability depending on the prevailing bacterial environment²⁶,²⁷. Different environmental conditions could induce distinct surface receptors phenotypes at the surface of bacteria, potentially changing dramatically their susceptibility to certain phages, resulting in entirely disparate outcomes in their interactions with a particular phage.
## Supplementary Figures

<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 1. Interactions between the three laboratory-adapted strains MG1655, BE, and BL21 with the phages of the *Antonina Guelin* collection.**</figcaption></figure>

Interactions between MG1655 and the collection were performed at a single phage concentration (5.10⁸ PFU/mL). Consequently, it is only possible to state whether a phage does not perform a lytic interaction (MLC score = Ø), performs lytic interactions with countable lysis plaques (MLC score = 1), or full clearing (MLC score = 2) at phage concentration of 5.10⁸ PFU/mL. For BE and BL21, the interactions were tested at all 3 phage concentrations, so the MLC score ranges from Ø to 4.

The number of phages infecting these strains was: BE, *n* = 45 (46.9%), BL21, *n* = 48 (50.0%), MG1655, *n* = 34 (35.4%).
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Supplementary Figure 2.** The phage-bacteria interaction matrix is nested and weakly modular.</figcaption></figure>

A. Rows and columns of the interaction matrix are rearranged using the BINMATNEST algorithm²⁸ in order to exhibit the matrix nestedness. Nestedness indicates that specific phages tend to infect the same bacteria as more generalist ones. The p-value is computed by generating 1,000 random matrices with the same size and average connectance as the original interaction matrix. B. Rows and columns of the interaction matrix are rearranged using the BRIM algorithm²⁹ to exhibit the matrix modularity. As in A, the p-value for the BRIM score is computed based on the generation of 1,000 null matrices with similar properties (size and connectance) as the original interaction matrix.
<figure><img src="image_16.png" /><figcaption>Supplementary Figure 3. RBP-encoding genes are hotspots of intra-genus genomic variability.</figcaption></figure>

Whole genome pairwise alignment of phage genomes belonging to the same genus, on five example genera (*Vectrevirus*, *Dhillonvirus*, *Gamaleyavirus*, *Felixounavirus*, and *Tequatrovirus*)³⁰. The stroke between two subsequent open reading frames (ORF) at the same locus indicates the percentage of amino acid identity between them. ORFs encoding for RBP are highlighted in light blue.
<figure><img src="image_17.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_18.png" /></figure>

**Supplementary Figure 4. Receptor-binding proteins are major determinants of the phage host range.**

A. RBP sharing is sufficient to explain host range similarity among phages, even among different genera. Left: Genomic map of the three selected phages exhibiting similar host ranges³⁰. AL505_P3 and T145_P2 belong to the same genus (*Vectrevirus*), whereas MT1B1_P3 is a *Kayfunavirus*. Genes that are conserved (> 30% amino acid similarity) are colored according to their functional annotation. Right: Host range of the three phages. Strains that which encode in their genome an ABC-dependent capsule with serotype K1 are indicated (black). For the three phages, plaque assay raw images are provided on three example bacteria (NILS67, ECOR-36, and ECOR-24) at the three MOIphage concentrations. B. T145_P2, AL505_P3, and MT1B1_P3 share a beta-propeller containing RBP, also found in the K1F phage. Predicted structure of the monomer of the conserved RBP among AL505_P3, T145_P2, and MT1B1_P3¹⁸. The RBP contains an endosialidase catalytic beta-propeller domain (PF12217, light blue) containing a nested beta-barrel (residues 104-543, light green) upstream of a beta-prism/tailspike domain (residues 544-718, pink) and variable N- and C-terminal extensions (gray and red). The X-ray diffraction structure of the K1F phage tailspike (TSP) monomer is also shown as a reference (PDB: 3GVJ). The name of each structural domain is given in the example of K1F TSP, and its position in AL505_P3 RBP #2, T145_P2 RBP #2, and MT1B1_P3 RBP is indicated by arrows.
<figure><img src="image_19.png" /><figcaption>Supplementary Figure 5. PDP110_P1 is likely to be a *Klebsiella* capsule K55-targeting phage.</figcaption></figure>

A. PDP110_P1 only infects three strains in the *Picard* collection. Interactions between phage PDP110_P1 (Siphovirus/Lambdavirus) and the *Escherichia* isolates belonging to the *Picard* collection. The name of each isolate infected by PDP110_P1 is provided. Isolates that encode a K55-serotype *Klebsiella* capsule are indicated (black triangle). Phage PDP110_P1 infects only three bacteria in the full *Picard* collection, including its isolation strain (0.7% of lytic interactions for this phage). Strikingly, the three infected bacteria were scattered across phylogroups A and C and matched exactly the three strains encoding a *Klebsiella* capsule with serotype K55. B. PDP110_P1 tends to generate halos on the strains it infects. Raw plaque assay results are shown for the three isolates that are infected by PDP110_P1 (NILS06, PDP110, and IAI34). The interactions are tested at three phage concentrations (5.10⁸, 5.10⁷, 5.10⁶ PFU/mL from left to right). Raw plaque assay results showed that each infection event led to the apparition of halos surrounding lysis plaques, which is commonly mentioned as an indicator of surface polysaccharide degradation by enzymes expressed by the phage¹¹. This strongly suggests that phage PDP110_P1 is a K55-targeting phage.
<figure><img src="image_20.png" /><figcaption>B</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>D</figcaption></figure>

<figure><img src="image_23.png" /></figure>

**Supplementary Figure 6. Bacterial Genome size significantly correlates with phage virulence upon infection, but the number of antimicrobial resistance (AMR) genes encoded in the bacterial genome does not.**

An LMM is fitted to the number of lytic interactions undergone by each bacterial strain (A-B) or the average MLC upon infection (C-D), taking as a covariate either the bacterial genome size (A-C) or the number of antimicrobial resistance genes encoded in the bacterial genome (B-D). The LMM takes the inverse relatedness matrix computed from the bacterial core phylogenetic distance as a random effect³¹. Significance is assessed based on the pMCMC which gives the probability of the estimated coefficient being null under the estimated posterior distribution.
<figure><img src="image_24.png" /><figcaption>**Supplementary Figure 7.** Coverage of the *Picard* collection by the most infecting *k*-uplet of phages of the *Antonina Guelin* collection. The number of productive lysis events performed by the most covering *k*-uplet of phages is computed for diverse values of *k* from 1 to 25 (blue). The gray dashed line corresponds to the most covering triplet (*k* = 3) of phages, which was chosen as a naive baseline to compete with our phage cocktail recommender algorithm. This triplet performs 63% of productive lysis on the *Picard* collection.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Supplementary Figure 8.** A four-step algorithm to recommend phage cocktails targeting *Escherichia* strains. The recommender algorithm for phage cocktails is designed as a four-step pipeline. Each step aims at determining whether each phage in the *Guelin* collection should be added to the recommended cocktail regarding a specific criterion which is defined by the input features and the algorithm that is used. Phages are ranked according to their predicted probability of infection, and top-ranked phages are added to the recommended cocktail when their probability is higher than 0.5. The algorithm goes sequentially through the four steps and stops whenever *k*=3 phages are recommended. In the recommendation process, the pipeline is forbidden to recommend two phages from the same genus and the same isolation strain to ensure phage diversity within the cocktail (see methods).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Supplementary Figure 9. Recommended cocktails can be partitioned into generic and tailored cocktails.</figcaption></figure>

A. Distribution of the steps of the pipeline at which each phage is recommended in the pipeline.

We recommend phage cocktails composed of three phages for each of the 100 test strains. The recommender algorithm goes sequentially through four steps. This design choice causes the algorithms to be partitioned into two categories: “Generic” and “Tailored”. The distribution of the steps of the pipeline at which phages are recommended, depending on the cocktail category, is represented. ** $p = 3.5 \times 10^{-12}$** (two-sided Mann-Whitney test). B. Number of haplotypes to which each cocktail is recommended. The haplotype of a bacterial strain is defined by the combination of its Sequence Type (ST), O-antigen serotype, and H-antigen serotype. The number of haplotypes is a measure of the genetic diversity of the bacteria that are recommended for a specific cocktail. Cocktails are partitioned between “Generic” and “Tailored” cocktails.
References

1. Tenaillon, B., O., Skurnik, D., Picard. The population genetics of commensal Escherichia coli. " *Nat Rev Microbiol* **8**, 207–217 (2010).

2. Jolley, K. A., Bray, J. E. & Maiden, M. C. J. Open-access bacterial population genomics: BIGSdb software, the PubMLST.org website and their applications. *Wellcome Open Res.* **3**, 124 (2018).

3. Bessonov, K. et al. ECTyper: in silico Escherichia coli serotype and species prediction from raw and assembled whole-genome sequence data. *Microb. Genomics* **7**, 000728 (2021).

4. Denamur, D., E., Clermont, O., Bonacorsi, S., Gordon. The population genetics of pathogenic Escherichia coli. " *Nat Rev Microbiol* **19**, 37–54 (2021).

5. DebRoy, C. et al. Comparison of O-Antigen Gene Clusters of All O-Serogroups of Escherichia coli and Proposal for Adopting a New Nomenclature for O-Typing. *PloS One* **11**, e0147434 (2016).

6. Ratiner, Y. A., Salmenlinna, S., Eklund, M., Keskimäki, M. & Siitonen, A. Serology and genetics of the flagellar antigen of Escherichia coli O157:H7a,7c. *J. Clin. Microbiol.* **41**, 1033–1040 (2003).

7. Nobrega, F. L. et al. Targeting mechanisms of tailed bacteriophages. *Nat. Rev. Microbiol.* **16**, 760–773 (2018).

8. Maffei, E. et al. Systematic exploration of Escherichia coli phage-host interactions with the BASEL phage collection. *PLoS Biol.* **19**, e3001424 (2021).

9. Liu, B. et al. Structure and genetics of Escherichia coli O antigens. *FEMS Microbiol. Rev.* **44**, 655–683 (2020).

10. Scholl, D. & Merril, C. The genome of bacteriophage K1F, a T7-like phage that has acquired the ability to replicate on K1 strains of Escherichia coli. *J. Bacteriol.* **187**, 8499–8503 (2005).

11. Beamud, B. et al. Genetic determinants of host tropism in Klebsiella phages. *Cell Rep.* **42**, 112048 (2023).

12. Haudiquet, M., Buffet, A., Rendueles, O. & Rocha, E. P. C. Interplay between the cell envelope and mobile genetic elements shapes gene flow in populations of the nosocomial pathogen Klebsiella pneumoniae. *PLoS Biol.* **19**, e3001276 (2021).

13. Nanayakkara, B. S., O’Brien, C. L. & Gordon, D. M. Diversity and distribution of Klebsiella capsules in Escherichia coli. *Environ. Microbiol. Rep.* **11**, 107–117 (2019).

14. Lam, M. M. C., Wick, R. R., Judd, L. M., Holt, K. E. & Wyres, K. L. Kaptive 2.0: updated capsule and lipopolysaccharide locus typing for the Klebsiella pneumoniae species complex. *Microb. Genomics* **8**, 000800 (2022).

15. Tesson, F. et al. Systematic and quantitative view of the antiviral arsenal of prokaryotes. *Nat. Commun.* **13**, 2561 (2022).

16. Kauffman, K. M. et al. Resolving the structure of phage-bacteria interactions in the context of natural diversity. *Nat. Commun.* **13**, 372 (2022).

17. Eddy, S. R. Accelerated Profile HMM Searches. *PLOS Comput. Biol.* **7**, e1002195 (2011).

18. Lin, Z. et al. Evolutionary-scale prediction of atomic-level protein structure with a language model. *Science* **379**, 1123–1130 (2023).

19. Pas, C., Latka, A., Fieseler, L. & Briers, Y. Phage tailspike modularity and horizontal gene transfer reveals specificity towards E. coli O-antigen serogroups. *Virol. J.* **20**, 174 (2023).

20. van Kempen, M. et al. Fast and accurate protein structure search with Foldseek. *Nat. Biotechnol.* **1–4** (2023) doi:10.1038/s41587-023-01773-0.

21. Stummeyer, K., Dickmanns, A., Mühlenhoff, M., Gerardy-Schahn, R. & Ficner, R. Crystal structure of the polysialic acid-degrading endosialidase of bacteriophage K1F. *Nat. Struct. Mol. Biol.* **12**, 90–96 (2005).

22. Boeckaerts, D., Stock, M., De Baets, B. & Briers, Y. Identification of Phage Receptor-Binding Protein Sequences with Hidden Markov Models and an Extreme Gradient Boosting Classifier. *Viruses* **14**, 1329 (2022).

23. Sørensen, A. N., Woudstra, C., Sørensen, M. C. H. & Brøndsted, L. Subtypes of tail spike proteins predicts the host range of Ackermannviridae phages. *Comput. Struct. Biotechnol. J.* **19**, 4854–4867 (2021).
24. Touchon, M. et al. Phylogenetic background and habitat drive the genetic diversification of Escherichia coli. *PLoS Genet.* **16**, e1008866 (2020).

25. Whitfield, C., Wear, S. S. & Sande, C. Assembly of Bacterial Capsular Polysaccharides and Exopolysaccharides. *Annu. Rev. Microbiol.* **74**, 521–543 (2020).

26. Rousset, F. et al. Genome-wide CRISPR-dCas9 screens in E. coli identify essential genes and phage host factors. *PLoS Genet.* **14**, e1007749 (2018).

27. Sato, M. et al. Expression of outer membrane proteins in Escherichia coli growing at acid pH. *Appl. Environ. Microbiol.* **66**, 943–947 (2000).

28. Rodríguez-Gironés, M. A. & Santamaría, L. A new algorithm to calculate the nestedness temperature of presence–absence matrices. *J. Biogeogr.* **33**, 924–935 (2006).

29. Barber, M. J. Modularity and community detection in bipartite networks. *Phys. Rev. E* **76**, 066102 (2007).

30. Gilchrist, C. L. M. & Chooi, Y.-H. clinker & clustermap.js: automatic generation of gene cluster comparison figures. *Bioinforma. Oxf. Engl.* **37**, 2473–2475 (2021).

31. Hadfield, J. D. MCMC Methods for Multi-Response Generalized Linear Mixed Models: The **MCMCglmm** R Package. *J. Stat. Softw.* **33**, (2010).