# Ecological ubiquity and phylogeny drive nestedness in phages-bacteria networks and shape the bacterial defensome

Chloé Feltin, Sylvain Piry, Benoît Moury, Lola Chateau, Karine Berthier, Jonathan M. Jacobs, Jules Butchacas, Linda Fenske, Lillian Ebeling-Koning, Theo H. M. Smits, et al.

► To cite this version:

Chloé Feltin, Sylvain Piry, Benoît Moury, Lola Chateau, Karine Berthier, et al.. Ecological ubiquity and phylogeny drive nestedness in phages-bacteria networks and shape the bacterial defensome. PLoS Pathogens, 2025, 21 (12), pp.e1013428. 10.1371/journal.ppat.1013428. hal-05441904

HAL Id: hal-05441904

https://hal.inrae.fr/hal-05441904v1

Submitted on 12 Jan 2026

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>
RESEARCH ARTICLE

# Ecological ubiquity and phylogeny drive nestedness in phages-bacteria networks and shape the bacterial defensome

Chloé Feltin¹*, Sylvain Piry¹, Benoit Moury¹, Lola Chateau¹, Karine Berthier¹, Jonathan M. Jacobs²,³, Jules Butchacas²,³, Linda Fenske⁴, Lillian Ebeling-Koning²,³, Theo H. M. Smits⁵, Cindy E. Morris¹, Clara Torres-Barceló⁶*

1 Pathologie Végétale, INRAE, Montfavet, France, 2 Department of Plant Pathology, The Ohio State University, Columbus, United States of America, 3 Infectious Disease Institute, The Ohio State University, Columbus, United States of America, 4 Bioinformatics and Systems Biology, Justus-Liebig-Universitat Giessen, Giessen, Germany, 5 Environmental Genomics and Systems Biology Research Group, Institute of Natural Resource Sciences (IUNR), Zurich University of Applied Sciences: ZHAW Zurcher Hochschule fur Angewandte Wissenschaften, Wädenswil, Switzerland, 6 PHIM Plant Health Institute, INRAE, University of Montpellier, CIRAD, Institut Agro, IRD, Montpellier, France

* clara.torresbarcelo@inrae.fr (CTB); chloe.feltin@inrae.fr (CF)

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

## Abstract

Identifying the ecological and evolutionary factors that shape phage-bacterial interactions is key to understanding their dynamics in microbial communities. Yet, such interactions remain poorly characterised in plant agroecosystems. Here, we investigate the ecological determinants of the interaction between a highly diverse set of 23 phages isolated from diseased apricot trees and 44 bacterial strains from the *Pseudomonas syringae* species complex collected either from diseased apricot trees, healthy plants or non-agricultural environment. Based on their ecological origin, we expected phages to preferentially infect bacterial strains from the same ecological context, forming modular host-range patterns. Contrary to these expectations, we discovered a significantly nested structure, suggesting generalised infection dynamics rather than local adaptation, primarily driven by the broad ecological dynamics of this pathosystem. Analysis of the bacterial genomes showed that both the profiles of anti-phage defence systems and the distribution profiles of prophages are strongly shaped by bacterial phylogeny. Furthermore, while the number of defence systems showed limited correlation with the breadth of bacterial sensitivity to phages, prophage abundance exhibited a strong, non-linear link with phage virulence. Together, these findings provide an ecological and evolutionary perspective on phage-bacterium infection networks and new insights into a better understanding of the role of phages in agricultural ecosystems.

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

OPEN ACCESS

**Citation:** Feltin C, Piry S, Moury B, Chateau L, Berthier K, Jacobs JM, et al. (2025) Ecological ubiquity and phylogeny drive nestedness in phages-bacteria networks and shape the bacterial defensome. PLoS Pathog 21(12): e1013428. https://doi.org/10.1371/journal.ppat.1013428

**Editor:** Breck A. Duerkop, University of Colorado Anschutz Medical Campus, UNITED STATES OF AMERICA

**Received:** August 5, 2025

**Accepted:** December 23, 2025

**Published:** December 30, 2025

**Copyright:** © 2025 Feltin et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data availability statement:** The phage-bacterium infection network, based on average values, has been deposited in the Viral Host Range Database hosted by Institut Pasteur (https://viralhostrangedb.pasteur.cloud/)
data-source/287/). The raw optical density (OD) kinetics data and the processed datasets used for analysis (e.g., AUC, inhibition values, phage effect probability) have been deposited in Zenodo (DOI: 10.5281/zenodo.17735985). The 39 bacterial genomes sequenced in this study have been deposited under the BioProject accession number PRJNA1292126.

**Funding:** The work described in our manuscript was founded by the Plant Health and Environment (SPE) Department of the French National Research Institute for Agriculture Food and the Environment (INRAE) (to C-TB) and Avignon University (AU) (to C-TB). CF received PhD funding from both the SPE Department of INRAE and AU. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

Viruses that infect bacteria, known as phages, are part of microbial communities and influence the abundance, diversity, and traits of their hosts. In an agriculture-related context, they are commonly considered as potential biocontrol agents, but studying the bases of fundamental phage-bacterial interactions may help us better understand the plant microbiome and its applications. Many factors influence these interactions, and identifying which ones matter most remains a challenge. In our study, we investigated how phages from diseased plants interact with bacteria collected from diseased and healthy plants, as well as from surrounding environments. We expected phages to mainly infect bacteria from similar environments, but instead observed that they often infected bacteria regardless of their source. This suggests that phage activity in this system has few barriers, reflecting the wide ecological distribution of their bacterial hosts. We further investigated how bacteria defend against phages by identifying both defence systems and prophages within their genomes and using this information to explore their contribution to bacterial resistance or sensitivity to phages. Together, our findings offer new insights into how phage-bacterium relationships evolve and function in plant ecosystems.

## Introduction

Phages profoundly influence bacterial diversity and population dynamics [1]. The impact of phages on bacterial ecology and evolution has been extensively examined in ecosystems like oceans and the human gut [2–7], whereas research in agricultural systems remains limited [8]. A deeper understanding of the ecological roles, population dynamics, and interactions of phages with bacterial communities in plant-associated systems, could provide fundamental insights to guide the development of phage-based biocontrol strategies and the sustainable management of microbial communities in agriculture.

Phage interactions with bacteria depend on a multitude of factors. The phage must recognise and infect a host, integrate or replicate, while avoiding or blocking bacterial defence systems [9,10]. A range of approaches have been employed to unravel the determinants of phage-bacterial interactions, including ecological and evolutionary analyses of natural populations and experimental evolution [11,12]. Mechanistic studies have identified bacterial receptors [13,14], with more recent work focusing on defence systems associated with susceptibility to phage [14–17]. For example, Vibrio lentus strains with identical phage receptors display varying levels of susceptibility depending on their phage defence elements [18]. However, research on prophages and defence systems in phytopathogenic bacteria remains limited and has yet to connect these elements to host range data [19–25].

The complexities of phage-bacterial interactions can be visualised as ecological networks, which exhibit distinct structural patterns, commonly modularity or nestedness [26]. Nestedness, a structure in which phages exhibit a gradient of host range
breadth and bacteria varying levels of susceptibility, is commonly observed with limited genetic diversity of phages or bacteria [26–28]. In contrast, modularity, where phages and bacteria interact within modules but rarely across them, tends to emerge at broader evolutionary or ecological scales, such as when interactions involve hosts across species, or genera, or are examined over extended time periods [11,29]. Nonetheless, modules can also evolve rapidly and emerge within species (e.g., [12]), although such events may be rare. A recent study introduced statistical methods that incorporate quantitative host–parasite interaction data into ecological networks [30]. By considering interaction strength rather than simple infection presence or absence, these approaches provide a more nuanced view of infection patterns, host susceptibility, and phage infectivity, capturing dynamics overlooked by binary analyses.

The *Pseudomonas syringae* species complex is particularly interesting for studying phage-bacterial interactions in agriculture, due to its high genetic diversity and broad ecological distribution across cultivated plants, wild environments and aquatic habitats [31–33]. *P. syringae* strains responsible for bacterial canker of apricot trees (*Prunus armeniaca*) are particularly diverse, spanning multiple phylogroups (PG01, PG02, PG03, PG07) [34]. It is hypothesized that the pathogenic success of these apricot strains may be driven primarily by general ecological fitness rather than specific virulence factors. These lineages are not apricot-specific but thrive across varied environments, raising questions about the ecological contexts of their phage interactions, whether within the host or external reservoirs [34]. Exploring phage infectivity patterns on these strains can reveal how phylogenetic and ecological diversity shapes phage-bacterium dynamics.

Historically, phage host range profiling was used to differentiate *P. syringae* strains [35–39]. More recent studies have emphasised the isolation and characterisation of novel *P. syringae* phages and their potential as biocontrol agents [40,41]. Most phage host range assays were confined to the isolation host and the intended target, reflecting an applied focus. Others included a broader set of strains, though typically limited to agricultural isolates. For example, phages infecting *P. syringae* pv. *porri* exhibited strict specificity to this variant, having been tested mainly against *porri* strains and a few non-*porri* isolates, and *P. syringae* pv. *actinidiae* phages were largely pathovar-specific, with some extending activity to closely related lineages [42].

We focused on the causal agents of bacterial canker of apricot trees to investigate the structure and ecological drivers of phage–bacterium infection networks. We hypothesised that phages would be locally adapted to their bacterial hosts, leading to a modular infection network structure that could be explained by the presence of prophages and diverse bacterial defence systems. To test this, we addressed the following questions: i) Are phages isolated from diseased apricot trees specifically adapted to resident *P. syringae* compared to closely-related strains from other environmental reservoirs? ii) Do closely-related bacteria isolated from different environments differ in their adaptation to phages? iii) What bacterial genomic traits underlie the observed infection network structure? We conducted an integrative approach utilizing a quantitative host-range assay involving 23 genetically-diverse *P. syringae* phages and 44 *P. syringae* strains. This was coupled with whole-genome sequencing of all bacterial strains for detailed annotation of prophages and defence systems. This comprehensive methodology provides mechanistic insights into how ecological origin and bacterial genomic content shape phage–bacterial interactions.

# Results

## Host range characterisation of a diverse panel of phages on *P. syringae* bacteria

We characterised the phage–bacterial infection network by testing 23 phages against 44 strains from the *P. syringae* species complex (hereafter referred to as *P. syringae* strains) in liquid culture over 48 hours of bacterial growth monitoring. Unlike traditional plaque assays, which primarily assess phage replication on solid media, this quantitative approach measures the extent of bacterial growth suppression in liquid conditions, thereby capturing a complementary aspect of phage–bacterium interactions [9,43]. Each interaction was quantified using two complementary metrics: 1) a continuous variable, corresponding to the percentage of bacterial growth inhibition by phages (0–100%) (S1A Table), and 2) a binary variable,
<figure><img src="image_4.png" /></figure>

indicating whether the inhibition was statistically significant (1) or not (0), based on triplicate measurements (S1B Table). Because these values can be interpreted from either the phage or bacterial perspective, we provide clear terminology to distinguish them. From the bacterial perspective, these values represent susceptibility (percentage inhibition) and sensitivity (binary significance of inhibition). Conversely, from the phage perspective, these values reflect virulence (percentage inhibition) and infectivity (binary success of infection).

Overall, the tested bacterial strains showed continuous distributions of susceptibility and sensitivity values, reflecting the diversity of interactions of the panel with phages (Fig 1A). On average, bacterial susceptibility to phages was moderate, with a 24% reduction in growth, and strains were sensitive to 12.3 out of the 23 tested phages (53.5%) (S1A Fig), indicating a broader sensitivity than suggested by inhibition levels alone. The most susceptible bacterial strains were DG1a.7, GAW0089, and DG12.6 (all from PG02), with average inhibition levels of 46%, 50%, and 42%, respectively (S1A Table). The most broadly sensitive strains, DG1a.7, 7C, and CCE0067 (all from PG02), were significantly inhibited by 20, 19, and 19 phages, respectively (87%, 83%, and 83% of tested phages respectively) (S1B Table). In contrast, bacterial strains CC1583, DG8.15, and DGA9CS5-8 (from PG10, PG07 and PG10 respectively) showed the highest resistance, with average inhibition percentages of only 3%, 6%, and 6%, respectively (S1A Table). DG8.15 (from PG07) resisted 17 phages (74% of tested phages) and seven other strains (19B, 41D, CST0086, DGA9CS5-15, FMU-107, and USA0088 from PG01, PG02, PG02, PG10, PG07 and PG02) were resistant to 16 phages (70%) (S1B Table).

From the phage perspective, average bacterial growth inhibition across the strain panel was 23%, reflecting moderate virulence. However, each phage was able to infect an average of 22.4 out of 42 tested *P. syringae* strains (53.3%) (S1B Fig), indicating a relatively broad infectivity despite moderate inhibition levels. The most virulent phages were Ghual01 (65% inhibition), Arace01 (57%), and Draal03 (44%) (S1A Table). The most infecting phages were Nican01, Arace01, and Ghual01 infecting 39 (93% of tested bacteria), 37 (88%), and 34 (81%) strains, respectively (S1B Table). Moreover, Arace01 and Ghual01 were the only phages capable of infecting strains from other species in the *Pseudomonas* genus. Arace01 infected a *P. rhizosphaerae* strain at a low efficiency of plating (EOP) of 0.005, whereas Ghual01 infected a *P. graminis* strain at a moderate efficiency (EOP=0.017). In contrast, Orimi01, Cygsa01, and Pavpe01 displayed the lowest virulence, with an average bacterial growth inhibition of 0%, 3%, and 4% respectively (S1A Table), while Orimi01, Cygsa01, and Touem01 had the lowest infectivity, infecting only 2 (5% of tested bacteria), 5 (12%), and 11 (26%) strains, respectively (S1B Table).

The relationship between the breadth of infection and the strength of inhibition can reveal underlying evolutionary dynamics, including trade-offs in both phages and bacteria [44,45]. When considering only significant inhibition interactions, we detected a positive correlation between the phage virulence and the phage infectivity (Spearman rho=0.56; p=0.01; S9C Fig), indicating that phages with broader host ranges tend to exhibit higher average virulence, thus not supporting the presence of a trade-off. However, this interpretation should be viewed with caution, as some degree of metric dependency may persist. In contrast, no correlation was observed between bacterial susceptibility and bacterial sensitivity (Spearman rho=0.24; p=0.12; S9D Fig), reflecting greater variability among bacterial strains in their response to phage infections and suggesting the potential for trade-offs in bacterial defensive strategies.

These findings reveal substantial variability in phage virulence and infectivity, as well as bacterial susceptibility and sensitivity, with over than half (53%) of combinations resulting in significant growth inhibition. While some phages exhibit a broad host range, extending even beyond *P. syringae*, others are highly specialised, only infecting a few strains, reflecting a complex interplay between phage infectivity and bacterial resistance.

## Nested structure of the phage-bacterial interaction network and association with bacterial traits

The structure of the phage-bacterial infection network was analysed using multiple algorithms to assess its nestedness and modularity. Given the diversity of both phages and bacterial hosts, spanning 13 phage genera and four phylogroups of *P. syringae*, and the ecological origins of phages (diseased apricot trees) and bacteria (diseased apricot trees, healthy
plants and other environmental reservoirs), we hypothesised a modular structure type. However, the network representing bacterial growth inhibition percentage by phages displayed a strong nested structure (Fig 1B), as indicated by scores of 0.48 and 0.44, on a scale varying from 0 to 1, with the WINE and wNODF algorithms, respectively (Table 1). To determine the statistical significance of the nested pattern, the network was compared to null models, which are simulations differing in the constraints applied during randomization, ranging from unconstrained (S, N) to fully constrained (B). Nestedness was highly significant (p<0.01) against nearly all null models, except null model B, which was shown to be excessively conservative [30].

Modularity analyses differed depending on the algorithm (Table 1). With the *spinglass* algorithm, selected for its good performance in terms of false positive rate [30], the modularity score was low (0.01, on a scale varying from -0.5 to 1) and was not significant, whatever the null model. The four alternative modularity algorithms produced slightly higher scores (0.09<modularity<0.25) and showed significant modularity (p<0.01) with both null models C1 and R1. These four algorithms show lower performances than *spinglass* in terms of false positive rates, but combined with null models C1 and R1, they show generally acceptable performances in terms of false positive and false negative rates [30]. These results suggest globally a slight degree of modular organisation within the network. The strong and significant nestedness and the weak modularity were further supported by the analysis of the binary data representing the significance of bacterial growth inhibition (S2 Table).

To assess whether the nested structure of phage-bacterial network is shaped by bacterial evolutionary history, we first tested for a phylogenetic signal in bacterial sensitivity summed by bacteria Fig 1. The analysis revealed a strong phylogenetic signal (λ=1.0; LR (λ=0) = 48.26; p<0.01), indicating that closely-related bacterial strains tend to be infected by a similar number of phages. To explore this pattern further, we evaluated the effect of phylogroup identity using both a Generalised Linear Model (GLM) with the bacterial sensitivity summed by bacteria (Fig 1C) and a PERmutational Multivariate ANalysis Of VAriance (PERMANOVA) on bacterial susceptibility profiles (S2A Fig). Both analyses confirmed a significant influence of phylogroup on bacterial susceptibility, explaining almost 10% of the variance (GLM: p<0.01; PERMANOVA: R²=0.09; p<0.01). In contrast, ecological origin had no significant effect in either test (GLM: p=0.87; PERMANOVA: R²=0.05; p=0.29). Specifically, PG02 strains emerged as the most sensitive to phages, being infected on average by 14.6 phages and showing the highest susceptibility, with 20.7% growth inhibition. PG01 strains showed intermediate sensitivity being infected by 10.4 phages and an average of 14.6% inhibition, followed closely by PG10 strains with 15.3% inhibition. In contrast, PG07 strains were the most resistant overall, infected on average by only 8.4 phages and displaying the lowest susceptibility with 8.73% inhibition.

The broad taxonomic diversity of the phages studied here, covering 13 genera, includes both specialists and generalists. To characterise phage host range, we calculated the Phylogenetic Host Range Index (PHRI), which incorporates

<table><caption>**Table 1.** Analysis of nestedness and modularity of the bacterial growth inhibition network combining 44 strains of the *P. syringae* species complex and 23 phages. Null models are described in Moury *et al.* [30]. For each null model and each algorithm, the frequency of simulations, among 100, showing a strictly higher nestedness (or modularity) score than that of the actual matrix is indicated. Significantly nested or modular networks (p≤0.05) are indicated in red.</caption><thead><tr><th rowspan="2">Analysis</th><th rowspan="2">Algorithm</th><th rowspan="2">Score (max.=1.0)</th><th colspan="8">Null model</th></tr><tr><th>B</th><th>N</th><th>C1</th><th>R1</th><th>S</th><th>C2</th><th>R2</th></tr></thead><tbody><tr><td rowspan="2">Nestedness</td><td>wNODF</td><td>0.44</td><td>0.44</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>WINE</td><td>0.48</td><td>1</td><td>0</td><td>0.01</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td rowspan="5">Modularity</td><td>spinglass</td><td>0.01</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0.07</td><td>0.29</td><td>0.18</td></tr><tr><td>edge betweenness</td><td>0.09</td><td>0</td><td>0.02</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td></tr><tr><td>fast greedy</td><td>0.24</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0.81</td><td>0.02</td></tr><tr><td>leading eigenvector</td><td>0.14</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0.94</td></tr><tr><td>beckett</td><td>0.25</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0.01</td><td>1</td></tr></tbody></table>

https://doi.org/10.1371/journal.ppat.1013428.t001
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Fig 1.** Networks of percentages of growth inhibition of 44 *P. syringae* species complex strains (rows) by 23 phages (columns). (A) Taxonomic organisation of the network presented on a maximum likelihood phylogenetic tree based on bacterial core-genomes and phage genera. A node with a maximum likelihood support value below 0.90 is marked with a red dot. (B) The same network after row and column permutations showing its nested pattern, with phylogroups and ecological origins of *P. syringae* strains displayed on the right side of panels A and B. (C) Bacterial sensitivity to phages by bacterial phylogroup and ecological origins. Letters indicate Tukey post-hoc test results from a Generalised Linear Model (GLM) with a Poisson distribution; ns indicates non-significant differences (p>0.05). (D) Phylogenetic Host Range Index (PHRI) across phage genera.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1013428.g001

the number, the evenness, and the phylogenetic diversity of host bacteria [46]. The PHRI values ranged from 0 (Orimi01, a specialist phage) to 0.38 (Arace01, a generalist phage), forming a gradient that reflects the nested structure of the phage-bacteria interaction network (Fig 1D). Based on PHRI classification, 9% of phages were categorised as specialists, 65% as generalists, and 26% as having an intermediate host range (Fig 1D). At the genetic level, significant differences in virulence were observed between phage genera (GLMM; p<0.01). The generalist genus *Aravirus* exhibited the highest virulence (52.8%), while the specialist genera *Cygnusvirus* and *Orionvirus* displayed markedly lower virulence levels, averaging 1.1% and 4.2%, respectively (S2B Fig). Additionally, phages with similar host infection profiles clustered by taxonomy, with those sharing >70% nucleotide identity (same genus) explaining 63.5% of the variation in virulence (PERMANOVA: R²=0.64, p<0.01), indicating genus-level host specificity (Fig 1A). For example, the two phages from the *Lyravirus sulafat* species and the *Ulliginvirus* genus exhibited highly similar infection host profiles.
<figure><img src="image_7.png" /></figure>

These results demonstrate that the phage-bacterial interaction networks are predominantly nested within the inves-
tigated panel of strains, with minimal modularity suggesting broad host range overlap among phages rather than local
adaptation based on ecological origin. Bacterial susceptibility is phylogenetically constrained, with PG02 strains being
highly sensitive and PG07 strains resistant. Additionally, although more than half of the phages were generalists, genus-
level similarity drives host-range patterns.

**Defence systems: profile and association with ecology and phylogeny**

The 44 strains of *P. syringae* analysed in this study were examined for their accessory genome content, with a partic-
ular focus on their repertoire of defence systems to counteract mobile genetic elements such as phages. Across these
strains, a total of 48 defence system types and 86 subtypes were identified (Fig 2A). The most abundant defence sub-
types included Septu (n=81), type IV Restriction-Modification (RM) (n=37), and Gabija (n=31) (Fig 2A). In contrast, 31
subtypes were sporadically detected across the bacterial genomes, including rare defence systems such as three Abi
subtypes (AbiH, AbiJ, and AbiO), two Asv subtypes (Asvl and Asvll), and three DTR subtypes (DTR6, DTR8, and DTR3)
(Fig 2A). Defence systems are present at an average rate of 9.93 per strain, with two strains (AF0015 and P3.01.09C09,
both from PG01) harbouring the highest number (15 systems), and four strains (DG1A_7 from PG02, CSZ0293, SZ0119,
and SZB0008 from PG10) the lowest number (5 systems) (Fig 2A). Each defence system can be classified into different
functional categories to better characterise its role. Specifically, 38.4% of defence system subtypes are associated with
abortive infection, 24.4% with nucleic acid degradation, and 37.2% remain unclassified (Fig 2A).

To assess whether the distribution of defence systems across bacterial genomes was influenced by phylogeny, we tested
for a phylogenetic signal in the total number of encoded systems per genome. A strong signal indicated that closely related
strains tend to harbour similar numbers of defence systems (λ=1.0; LR test (λ=0): 44.18; p<0.01) (S3 Fig), as described in a
previous study [24]. Phylogroup significantly explained variation in defence system profiles (PERMANOVA, R²=0.14, p<0.01),
with a weaker effect from ecological origin (R²=0.06, p<0.01), while growth rate of the bacteria had no significant impact
(R²=0.02, p=0.51) (Fig 2A). Post-hoc comparisons based on this analysis revealed that strains from PG10 displayed signifi-
cantly distinct defence system profiles compared to those from PG02 and PG07, while strains isolated from healthy plants
exhibited distinct defence system profiles compared to other ecological origins. However, none of these factors significantly
affected the total number of defence systems, although the effect of phylogroup approached significance (GLM: p=0.06)
(Fig 2B). This phylogroup trend is supported by a larger dataset of 590 P. syringae genome assemblies, where PG01 strains
encode the highest number of defence systems, while PG07 and PG10 harbour some of the lowest counts [24].

While closely-related strains tend to possess similar numbers of defence systems, we further investigated whether they
also share the same subtypes of defence systems. To this end, we calculated the percentage of shared defence system
subtypes between each pair of bacterial strains, where 100% indicates complete identity in the sets of subtypes. A weak
but significant positive correlation was observed between the percentage of shared defence system subtypes and the
genome-wide Average Nucleotide Identity (ANI) among bacterial strains (GLMM; estimate = 0.015; p < 0.01) (S4 Fig). This
relationship suggests that some defence systems are partly inherited through vertical transmission, reflecting the influence
of phylogeny on their distribution. Furthermore, only 1.14% of the identified defence systems were linked to prophages
within bacterial genomes, indicating that the vast majority are located outside prophage regions. This low prevalence
indicates that while defence systems may spread horizontally among bacterial populations, their dissemination is likely
mediated by mobile genetic elements other than prophages.

Defence systems may impose varying energetic costs on bacterial cells, potentially influencing their allocation of
resources [47,48]. In strains of *P. syringae*, defence systems related to nucleic acid degradation comprised, on aver-
age, 54% of the total defence systems, while abortive infection systems, which are considered more metabolically costly
[49,50], constituted an average of 25% (S3 Fig). We then revealed that the proportion of defence system categories were
significantly shaped by phylogroup, explaining 16% of the variance (PERMANOVA; R² = 0.16; p = 0.03), and growth rate
<figure><img src="image_8.png" /><figcaption>**Fig 2.** Defence system distribution in bacterial genomes and relationship with bacterial sensitivity to phages. (A) Matrix showing the number of defence systems present in the 44 *P. syringae* strains, organised on a maximum likelihood phylogenetic tree based on their core-genomes. A node with a maximum likelihood support value below 0.90 is marked with a red dot. The phylogroup affiliation of each bacterial strain is indicated using colours, and defence systems are represented by subtypes with colours indicating their respective functional category. (B) Number of defence systems per bacterial strain according to their phylogroups and ecological origins. ns: not significant (p > 0.05), according to post-hoc Tukey tests from a Generalised Linear Model (GLM) with a Poisson distribution. (C) Bacterial sensitivity to phages as a function of the number of defence systems, with linear regressions performed within each phylogroup.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1013428.g002

(R²=0.07; p=0.04), but not ecological origin (R²=0.08; p=0.11). Pairwise comparisons based on this analysis showed that strains from PG07 had significantly different defence category profiles compared to those from PG01 and PG02.

Our results revealed that the defence arsenal profiles in *P. syringae* are shaped by phylogeny and slightly by ecological origin, consistent with partial vertical inheritance of defence system profiles. Overall, these findings underscore the complexity of bacterial defence strategies and indicate that multiple evolutionary and ecological factors influence phage resistance.

## Prophages: profile and association with ecology and phylogeny

We identified 88 prophages among the 44 genomes of *P. syringae*. On average, each strain carried two prophages, with one strain (FMU-107) harbouring five prophages and 18 strains (41%) containing a single cryptic prophage (Fig 3A).
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig 3.** Diversity and distribution of prophages in bacterial genomes and relationship with bacterial sensitivity to phages. (A) Distribution of prophages in the genomes of the 44 *P. syringae* strains, categorised by genus and family-like cluster (VipTree). Bacterial strains are organised based on a maximum likelihood phylogenetic tree based on their core-genomes, with their phylogroups and ecological origins indicated. A node with a maximum likelihood support value below 0.90 is marked with a red dot. (B) Abundance of prophages per bacterial strain according to their phylogroups and ecological origins. ns: not significant (p>0.05), according to post-hoc Tukey tests from a Generalised Linear Model with a Poisson distribution. (C) Bacterial sensitivity to phages as a function of the number of prophages, coloured by phylogroups.</figcaption></figure>

https://doi.org/10.1371/journal.ppat.1013428.g003

Among these prophages, 36% were classified as high quality (complete or nearly complete according to CheckV software), 39% as medium quality, and 25% as low quality. Prophage sizes ranged from 12.0 kb to 89.8 kb, with an average size of 32.2 kb. The size distribution exhibited three peaks at approximately 18 kb, 44 kb, and 57 kb. To classify these
<figure><img src="image_11.png" /></figure>

prophages taxonomically, we used VIRIDIC, which identified 65 novel phage species and 42 novel genera (S3 Table), while proteomic tree analysis with Viptree revealed six distinct prophage clusters, corresponding approximately to family-level taxonomic divisions (Fig 3A) [51]. One of these, the CL4 family-like prophage, was found in all 44 *P. syringae* strains, which span four distinct phylogroups (Fig 3A). CL4 family-like prophage is a tailocin, specifically an R-type syringacin, a bacteriocin derived from a phage [52,53]. These tailocins are widely distributed across the *P. syringae* species complex and are thought to play a key role in inter-bacterial competition by enabling strains to kill closely-related competitors [54,55]. Moreover, we found that the taxonomy of the CL4 family-like prophage mirrors the phylogeny of their bacterial hosts across 68 strains from 10 phylogroups (including the 44 strains analysed in this study), with a significant cophylogenetic congruence signal (Procrustes Approach to Cophylogeny (PACo): residual sum of squares = 0.15; *p* < 0.01) (S5 Fig).

The abundance of prophages (i.e., the number of prophages per bacterial genome) and the prophage distribution profiles in bacterial genomes are strongly shaped by phylogeny. A significant phylogenetic signal was detected in prophage abundance (λ = 1.0; LR (λ = 0) = 46.48; *p* < 0.01), indicating that closely related strains tend to harbour similar numbers of prophages (Fig 3B). Phylogroup significantly explained 16% of variation in prophage distribution profiles (PERMANOVA, R² = 0.16; *p* < 0.01), with a weaker effect from ecological origin (R² = 0.07; *p* < 0.01), while growth rate had no significant impact (R² = 0.03; *p* = 0.09) (Fig 3A). All distinct phylogroups, as well as ecological origins, exhibited significantly different prophage distribution profiles, supporting the idea that lineage and ecological context influence the genomic landscape of prophages. However, none of these factors significantly affected the abundance of prophages (GLM: phylogroup, *p* = 0.27; ecological origin, *p* = 0.92; growth rate, *p* = 0.30), despite some observed variation. Indeed, PG07 strains showed the highest prophage abundance, nearly double that of PG10 strains (Fig 3B). However, this pattern contrasts with findings from broader dataset (590 *P. syringae* genome assemblies), in which PG02 strains encoded among the highest numbers of prophages, PG07 among the lowest, and PG01 and PG10 showed intermediate levels [24]. In contrast, differences across ecological origins were minimal (Fig 3B).

In conclusion, the comprehensive analysis of prophages within the *P. syringae* species complex revealed a diverse array of complete and incomplete prophages integrated into the genomes. The prophage distribution profiles are strongly explained by bacterial phylogeny and ecological origins of the host strains.

## Contribution of bacterial defence systems and prophages to bacterial sensitivity to phages

Since bacterial defence systems target phages [50], a negative correlation is expected between their number and bacterial sensitivity to phages. To test this, we examined the correlation between bacterial sensitivity and the number of encoded defence systems, while accounting for phylogenetic background. When considering all bacterial strains together, no overall correlation was identified (*p* = 0.36) with a phylogroup effect (*p* < 0.01) and no effect of ecological origins (*p* = 0.88). For both PG01 and PG02, no significant correlation was detected (linear regression for PG01: adjusted R² = 0.12; *p* = 0.18; for PG02: adjusted R² = -0.07; *p* = 0.94) (Fig 2C). In contrast, a negative correlation was detected in PG10, supporting our expectation (linear regression for PG10: adjusted R² = 0.34; *p* = 0.05) (Fig 2C). In contrast, PG07 exhibited an unexpected positive correlation, suggesting that strains with a higher number of defence systems were more sensitive to phages (linear regression for PG07: adjusted R² = 0.92; *p* < 0.01) (Fig 2C).

Although prophages can contribute to bacterial immunity against phage infection [56], we found no significant association between prophage abundance and bacterial sensitivity to phages (*p* = 0.82), and neither phylogroup (*p* = 0.31) nor ecological origin (*p* = 0.80) had a detectable effect (Fig 3C). Moreover, no correlation was observed between prophage abundance and the number of defence systems encoded in bacterial genomes, suggesting that these systems do not provide protection against prophage integration or maintenance (linear regression; adjusted R² = -0.02; *p* = 0.64).
Influence of multiple bacterial and phage traits on bacterial growth inhibition

Having identified several bacterial and phage factors that individually influenced the bacterial growth inhibition, we next sought to integrate these explanatory variables into a unified framework to gain a more comprehensive understanding of phage-bacteria interaction dynamics. The primary model included prophage abundance and the total number of defence systems per strain as bacterial traits associated with phage resistance, alongside phage- and bacterium-related traits examined in this study. We also tested an alternative model incorporating PCA-derived variables related to bacterial defensome, but it did not improve model performance (Akaike Information Criterion (AIC) = 16.70) compared to the initial model (AIC = 16.69). Therefore, we retained the first model for its greater interpretability.

Consistent with the above findings, bacterial phylogeny and phage taxonomy emerged as the strongest explanatory variables for inhibition percentage, with an importance weight of 1.0 across the 256 models (Table 2). In addition, a variable indicating whether the phage had been isolated using one of the tested bacterial strains was significantly associated with inhibition level. Specifically, phages isolated from strains within the tested host range exhibited a slightly lower inhibition capacity (mean = 15.1%) compared to those isolated from external strains (mean = 17.7%).

The model identified the total number of prophages in bacterial genomes as a significant predictor of inhibition of bacterial growth (sum of weights = 0.72; Table 2). Inhibition levels varied nonlinearly with prophage abundance: inhibition levels were 15.6%, 14.7%, 17.5%, 28.4%, and 10.9% for strains carrying 1–5 prophages, respectively. While prophage abundance did not predict bacterial sensitivity, it appeared to influence bacterial susceptibility, suggesting that prophage-rich bacteria may be more susceptible to highly-effective phages. In contrast, the total number of defence systems did not significantly explain bacterial growth inhibition (sum of weights = 0.38; Table 2), despite its correlation with the bacterial sensitivity for some phylogroups.

Overall, these findings highlight that the primary determinants of the bacterial growth inhibition by phage infection within the *P. syringae* species complex are bacterial phylogeny and phage taxonomy rather than bacterial ecological origins. The integration of multiple explanatory variables within a unified modelling framework further revealed that prophage abundance, but not the total number of defence systems, significantly contributes to variations in inhibition levels. This suggests that while prophages may not act as strict barriers to infection, their presence can modulate the efficacy of phage infection. Together, these results emphasize the complex and multifactorial nature of bacterial susceptibility to phages, shaped by both evolutionary history and prophage-mediated genomic landscapes.

<table><caption>Table 2. Relative importance of explanatory variables influencing bacterial growth inhibition by phages. Sum of AIC weights was computed from model averaging across the top-ranked models whose cumulative weights reached ≥ 0.95, selected from a full set of 256 candidate models generated by multi-model inference. The relative importance of each variable corresponds to the sum of AIC weights across all retained models in which it appears. Bold values indicate variables with high relative importance (sum of weights ≥ 0.60), selected as significant predictors in the analysis. N Models: Number of retained models in which each variable appears. Variables related to phages are marked with the symbol φ; bacterial variables are unmarked.</caption><thead><tr><td>Explanatory variable</td><td>Sum of AIC weights</td><td>N Models</td></tr></thead><tbody><tr><td>Phage genera φ</td><td>1.00</td><td>54</td></tr><tr><td>Bacterial phylogroups</td><td>1.00</td><td>53</td></tr><tr><td>Number of prophages in bacterial genomes</td><td>0.72</td><td>29</td></tr><tr><td>Phages isolated with bacteria from the data set φ</td><td>0.67</td><td>31</td></tr><tr><td>Number of defence systems in bacterial genomes</td><td>0.38</td><td>26</td></tr><tr><td>Bacteria from the data set having isolated phages</td><td>0.37</td><td>24</td></tr><tr><td>Bacterial growth rate</td><td>0.29</td><td>23</td></tr><tr><td>Bacterial ecological origins</td><td>0.12</td><td>21</td></tr></tbody></table>

https://doi.org/10.1371/journal.ppat.1013428.t002
<figure><img src="image_12.png" /></figure>

Discussion

Our study uncovered a strong nested structure in the phage-bacterial infection network involving genetically-diverse *P. syringae* strains and phages. This discovery is in strong contrast with the modular structure frequently reported [11,29]. This pronounced nestedness likely reflects the ubiquitous nature of the highly diverse *P. syringae* complex, mirrored by a phage community ranging from specialists to generalists in host range. While bacterial phylogeny had a marked influence on the nestedness of the infection network, we found no evidence of local phages-bacteria co-adaptation. Defence systems and prophages contributed to phage resistance, but their distribution was largely shaped by bacterial phylogeny and ecological origin, rather than phage pressure alone. Together, these findings highlight how evolutionary and ecological forces shape phage-bacterial infection networks in natural environments.

Ecological origin and phylogenetic influences on bacterial sensitivity to phages

The observed nested infection network is unexpected when we consider the high genetic diversity of the phages (13 genera, 17 species) and bacteria (4 phylogroups corresponding to species levels and 7 clades), as well as the varied ecological origins of the bacterial strains included in this study. Such a diversity of organisms could be expected to reveal local adaptation, where phages perform better on sympatric than allopatric bacterial hosts, in line with what has been demonstrated previously [11,57]. However, our findings suggest that phage-host interactions in this pathosystem transcend ecological boundaries. Phages isolated from apricot tree soil infected bacterial strains from both diseased orchards and aquatic environments with comparable efficiency, and bacterial origin did not broadly influence bacterial sensitivity/susceptibility to phages. This absence of modularity, combined with a lack of ecological signal in host range, points to broader ecological connectivity and population flow within the *P. syringae* species complex. These results align with previous work [58] showing that *P. syringae* population structures are shaped more by fluctuating ecological dynamics than by geographic isolation, likely due to documented exchanges between agricultural and aquatic environments [59]. Overall, the nested structure observed here implies that interaction dynamics in this system are maintained by widespread generalism and ongoing ecological mixing, rather than niche-specific adaptation.

Bacterial phylogeny plays a central role in shaping both the composition of defence systems, the distribution of prophages, and the sensitivity of *P. syringae* strains to phages. In our study, PG02 strains, known for their ubiquitous presence across habitats and for containing strains capable of infecting a broad range of plant species [33,60], were also the most sensitive to phages, carried a moderate number of defence systems, and encoded one of the highest prophage loads [24]. In contrast, PG07 strains, which are predominantly associated with diverse environmental reservoirs [60], were the most resistant to phages despite encoding relatively few defence systems and prophages. PG01 strains, mainly isolated from diseased plants [60], showed intermediate bacterial sensitivity to phages, possessed the greatest number of defence systems, and encoded a moderate prophage load [24]. Finally, PG10 strains, exclusively environmental, showed moderate bacterial sensitivity to phages, low numbers of defence systems, and a moderate prophage count [24]. These findings suggest a potential trade-off between ecological ubiquity and bacterial sensitivity to phages, a hypothesis that requires direct experimental investigation. In a similar vein, a recent study reported that effector gene counts in *P. syringae*, which are key to pathogenicity, positively correlate with phage defences and show partial conservation within phylogroups [24].

Phylogroups show consistent patterns in defence system profiles, with closely-related strains sharing portions of their
defence arsenals, while ecological origins accounted for only minor similarities. Strains within the same phylogroup also
tend to retain consistent features in their defensome, including the number and functional categories of defence systems.
Collectively, these findings suggest that each phylogroup may adopt a common anti-phage strategy, implemented through
distinct combinations of defence systems, mainly acquired by horizontal gene transfers [50,61]. For strains belonging to
phylogroup PG10, we observed a negative correlation between the number of defence systems and sensitivity to phages,
<figure><img src="image_13.png" /></figure>

suggesting a protective effect from accumulating diverse systems in their genome [17]. Conversely, strains in PG07, although generally the most resistant to phages, showed a positive correlation between bacterial sensitivity to phages and defence systems. This counter-intuitive pattern could reflect historical predation by phages, where strains subjected to intense phage attack would have accumulated defence systems [62]. Alternatively, the apparent lack of protection may arise because these systems target phages or mobile genetic elements (MGEs) other than those tested here. However, it is important to interpret these trends cautiously, as they are based on relatively small sample sizes (n=6 for PG07; n=10 for PG10), and our understanding remains partial given the continual discovery of novel defence systems [63]. Although the sensitivity of specific phylogroups to phages may be correlated to their defence systems, the overall contribution of these systems to phage resistance in *P. syringae* strains appears limited and may depend on the activity of particular defence mechanisms [64].

## Limitations and alternative approaches for understanding phage-bacterium interactions

Liquid culture assays offer the advantage of measuring a proxy for phage virulence by continuously monitoring bacterial growth kinetics in a high-throughput format. This automated approach enhances reproducibility, allows the simultaneous processing of a large number of samples, and enables the use of diverse metrics to quantify phage-induced inhibition of bacterial growth [65], including host-range, and potential indicators of bacterial resistance to phages [9,43]. In our dataset, certain phage-bacterium combinations exhibited signs of evolutionary rescue (when a bacterial population initially suppressed by phages recovers through the growth of resistant mutants) [66] (S11 Fig). Future studies could build on these findings by examining their genetic co-evolution through temporal sampling and genomic analyses. Nevertheless, we acknowledge that the phenomenon of lysis from without (referring to bacterial lysis without phage replication, typically resulting from high phage adsorption rates or abortive infections) cannot be entirely ruled out in our dataset [67]. Follow-up experiments measuring phage titres, bacterial population dynamics, or adsorption efficiency would be valuable to confirm the underlying mechanisms of bacterial suppression. Such analyses would also refine our understanding of the role of defence systems during phage infection [14,17]. In particular, quantifying phage adsorption rates could help discriminate between pre- and post-adsorption resistance: successful adsorption without subsequent lysis would suggest the involvement of intracellular defence mechanisms, whereas failed adsorption would point to surface receptor incompatibility as the main barrier to infection.

Analysis of bacterial membrane receptors recognised by phages is therefore essential in interpreting host range patterns [12,68]. Such factors have enabled accurate prediction of phage-bacteria interactions at the strain level in *E. coli* [14]. However, receptors for phages infecting *P. syringae* remain poorly characterised compared to other bacterial models, with only three distinct receptors identified in the *P. syringae* phage receptor database compared to over 50 in *E. coli* [69]. Increasing research on phage-phytopathogen interactions will help to close this knowledge gap and improve our understanding of the molecular basis underlying the observed interaction patterns.

A complementary strategy to uncover bacterial genetic determinants of sensitivity to phages and resistance phenotypes involves genotype-phenotype association analyses, including approaches based on gene presence-absence matrices, phylogenetic trait correlation methods and genome-wide association studies (GWAS). This latter method has recently been applied to 80 *P. aeruginosa* strains, identifying host genes that influence susceptibility to two phages, particularly those related to flagella and LPS/O-antigens associated with lipopolysaccharide biosynthesis [70]. While our dataset of 44 *P. syringae* genomes may provide a sufficient basis for GWAS, challenges remain. The high genomic diversity of the *P. syringae* species complex and the large number of phages tested (n=23), each potentially interacting with different receptors or resistance pathways, add considerable complexity. Future studies may address this by focusing on specific phage-host combinations or by applying comparative genomics approaches that can help disentangle the genetic basis of resistance across diverse strains.
<figure><img src="image_14.png" /></figure>

## Applied perspectives and future directions

From an applied perspective, we identified several generalist phages capable of infecting a wide range of *P. syringae* species complex strains. Understanding the molecular basis of their generalism could provide important phage evolutionary clues, e.g., whether generalism involves flexible receptor-binding proteins, the ability to evade diverse bacterial defence systems, or other mechanisms. Many *P. syringae* strains are non-pathogenic members of the plant microbiome, highlighting the need for targeted phage cocktails that selectively control pathogenic strains while preserving beneficial microbes. Such specificity offers a more sustainable and ecologically-balanced approach to biocontrol compared to using broad-spectrum generalist phages [71]. This exemplifies how fundamental ecological understanding of phage-bacterium infection networks can directly inform and enhance the development of effective biocontrol strategies.

To further test our hypothesis that ecological connectivity shapes the nested structure of phage-P. syringae interactions, future studies should isolate phages from diverse non-agricultural environmental, such as water and healthy plants to include cases of true sympatry. These newly isolated phages should then be tested against the same bacterial strains used in our study, or against new strains from corresponding habitats. If similar nested network structures emerge, it would strengthen the idea that phage and *P. syringae* populations continuously migrate between agricultural and non-agricultural environments, forming an ecologically-connected meta-community. Moving forward, several key research priorities should guide future investigations: verification of the role of specific defence systems through genetic manipulation like gene deletion and heterologous expression; systematic investigation of receptor usage by our studied phages to reveal molecular determinants of host range; and strategic expansion of the bacterial strain collection to improve phylogenetic and ecological coverage, with particular focus on understudied phylogroups. These complementary approaches would collectively strengthen our understanding of the molecular and ecological factors driving the phages-bacteria interactions observed in this system. Extending such analyses to other ecologically and phylogenetically-distinct plant-associated bacteria, such as *Xanthomonas*, *Ralstonia*, or *Agrobacterium*, could reveal whether the patterns observed in *P. syringae*-phages infection network are generalisable or reflect unique evolutionary dynamics of this highly diverse species complex.

## Materials and methods

### Phages panel

The 23 phages infecting strains of the *P. syringae* species complex used in this study were isolated from soil samples collected in apricot orchards affected by bacterial canker and fully described in a previous study [72]. Briefly, this phage collection exhibits high genetic diversity, encompassing 13 genera and 17 species, representing the three most common morphological types (myovirus, podovirus, and siphovirus) and comprising 20 virulent and 3 temperate phages. The collection included two characterised DNA packaging mechanisms (direct terminal repeat (DTR) and headful) alongside two newly identified but unverified mechanisms, while the remaining packaging strategies remain unknown. The genomic GC content of the phages ranged from 46% to 58%. All phage stocks were standardised to a concentration of 10⁷ PFU/mL and stored in SM buffer (100 mM NaCl, 8 mM MgSO₄·7H₂O, 50 mM Tris-HCl, pH 7.4) at 4°C.

### Bacterial strains panel

The 44 *P. syringae* species complex strains were obtained from the Plant Pathology Unit collection at INRAE (Montfavet, France) and were not isolated as part of this project (S4 Table). The strains belonged to seven clades within four phylogroups (PG01, PG02, PG07, PG10) corresponding to the species levels, according to Berge *et al.* [60]. The strains were sampled from three different environments: diseased apricot trees (exhibiting bacterial canker symptoms) [34], healthy plant tissues, and non-agricultural environments (S4 Table). Bacterial strains were revived from -80°C storage and cultured on King's B agar [73]. Cultures were maintained at 4°C for a maximum of one week.
<figure><img src="image_15.png" /></figure>

## Sequencing and analysis of whole genomes of bacterial strains

The bacterial cultures were grown on King's B (KB) agar medium for 48 hours at 25°C. A bacterial suspension was then prepared to achieve an approximate concentration of 2 × 10⁹ CFU in 4 mL of demineralised water. The suspension was subsequently centrifuged at 3,215 g for 10 minutes at 20°C to pellet the bacterial cells, after which the supernatant was carefully discarded. Genomic DNA was extracted using the Qiagen Genomic-tips kit (Qiagen, Valencia, CA, USA) following the manufacturer's protocol.

Of the 44 *P. syringae* strains, 39 were sequenced using 20 µL of genomic DNA, with sequencing performed by Plasmidsaurus Inc (San Francisco, USA) using Oxford Nanopore Technology (ONT) and a custom analysis and annotation pipeline (https://plasmidsaurus.com/). The long-read sequencing library was prepared without amplification using v14 library preparation chemistry, and sequencing was conducted with R10.4.1 flow cells. The genomes of five additional strains were obtained from the RefSeq library at NCBI: CC1559 (GCF_000452685.1), USA0007 (GCF_000452545.1), 41A (GCF_000935775.1), MAFF302273<sup>PT</sup> (also known as JD08) (GCF_000145925.1), and CC1583 (GCF_000452665.1).

The core genome analysis was performed using a private database within the EDGAR 3.5 platform [74]. The core genome was defined as the set of orthologous genes present in all 332 genomes of the *P. syringae* species complex, including the 44 strains analysed in this study. A total of 1,346 core genes were identified per genome, corresponding to 446,872 aligned genes across all strains. Each set of orthologous genes was individually aligned at the protein level using MUSCLE [74]. The resulting alignment were then concatenated into a single alignment, yielding approximately 512,223 amino acid residues per genome and 170,058,036 in total across all genomes. Phylogenetic relationships were inferred using the FastTree algorithm within EDGAR, which constructs maximum-likelihood phylogenetic trees based on the core genome alignment [75]. All phylogenetic trees of bacterial strains shown in the figures were rooted at the midpoint of the longest path and the nodes with a maximum likelihood support value below 0.90 were marked with red dots.

## Identification of defence systems and prophages in bacterial genomes

The defence systems of the 44 *P. syringae* strains were identified using *DefenseFinder* version 1.2.2 [76] based on assembled nucleotide sequences from the bacterial genomes. Each defence system type was categorised into one of three functional groups: nucleic acid degradation, abortive infection (Abi), or undetermined (ND).

The prophages within the 44 *P. syringae* strains were identified using *geNomad* version 1.5.0 [77] based on assembled nucleotide sequences of bacterial genomes. As prophage prediction tools are still imperfect and may miss or misclassify certain elements, these results should be interpreted with caution. Among the 97 identified prophage sequences, only those classified as viral by *geNomad* and *CheckV* [78] were retained. Specifically, sequences labelled as "No Terminal Repeats" in the topology column and those with a "viral_genes" value of 0 were excluded. This filtering resulted in a final dataset of 88 verified viral prophages, accounting for 91% of the total identified prophages across the 44 strains. Additionally, a genetic analysis of the prophages was performed using *VIRIDIC* [79] to classify them into genera and species (S3 Table). *VipTree* was used to generate prophage clusters based on a proteomic tree analysis using distance-based methods (BIONJ from tBLASTx scores) [80]. *geNomad* and *VipTree* were then employed to identify prophages belonging to the CL4 cluster and to construct a proteomic tree encompassing 68 CL4 prophages (24 reference strains + 44 analysed strains) [81].

To assess the congruence between the taxonomy of CL4 prophages and the phylogeny of their bacterial hosts, we performed a PACo (Procrustes Approach to Cophylogeny) analysis. Phylogenetic distance matrices were generated from the host and prophage trees using the *cophenetic* function [82]. These matrices were then centred and corrected for negative eigenvalues using the "Cailliez" correction method. Congruence was tested using the *PACo* function from the *paco* R package with 1,000 permutations to assess statistical significance. The co-phylogeny visualisation was performed with *phytools* R package and *cophylo* function [83].
<figure><img src="image_16.png" /></figure>

## Phage host range analysis

**a. Experimental setup.** To analyse the host range of phages infecting *P. syringae*, bacterial growth kinetics in the presence of phages were measured using 96-well plates [65]. This liquid culture approach has been compared to other methods, including qPCR, spot assays, and plaque assays, demonstrating its relevance for estimating phage host range across diverse bacterial-phage combinations [65,84–86]. Given the broad panel of phages and bacterial strains analysed, we selected this method as it provides an optimal balance between throughput, reproducibility, and comparability, while avoiding the more labour-intensive phage progeny quantification assays.

The dataset comprised 44 bacterial strains tested in biological triplicate without phage as a positive control, and the same 44 strains in biological triplicate exposed individually to each of the 23 phages (S6 Fig). Each phage was tested against each biological triplicate of bacterial strain, resulting in 3,036 growth kinetic measurements (23 phages × 44 bacterial strains × 3 replicates) (S6 Fig). Negative controls (culture medium) were used to check for contamination, with plates containing more than two contaminated wells out of four discarded (S6 Fig and S7A Fig). Two internal controls were included on each plate: the *P. syringae* strain P1.01.01B03 (PG01a, isolated from apricot buds in France in 2014) and phage Ghuch01, chosen because it completely inhibited the growth of strain P1.01.01B03 under the experimental conditions used in this study [34,72] (S6 Fig and S7B-C Fig). A plate was considered valid only when all three controls were visually confirmed: the internal bacterial control exhibited visible turbidity, while both the internal phage control and the media controls showed clear wells (S6 Fig).

Bacteria were incubated in 5 mL of liquid KB medium at 25°C and 200 rpm for 4 hours. Each replicate was standardised to an $OD_{600nm} = 0.01$ (~$10^7$ CFU/mL) in phosphate buffer and diluted in liquid KB medium to a final concentration of $10^5$ CFU/mL. Each 96-well spectrophotometric plate (NUNC, Thermo Scientific 267334) contained most phages confronted with a single bacterial strain, with positive control wells containing 180 µL of bacterial suspension and 20 µL of SM buffer. Experimental wells comprised 180 µL of bacterial suspension and 20 µL of phage suspension at $10^7$ PFU/mL, corresponding to a multiplicity of infection (MOI) of 10, as determined from preliminary optimization tests. Positive control measurements included between two and seven technical replicates per biological replicate (except for bacterial strains 41A and 93D, which had only one technical replicate per biological replicate). Measurements were taken using the SPECTROstar Omega spectrophotometer (BMG Labtech, France) every 5 minutes for 48 hours at 25°C, with intermittent orbital shaking at 200 rpm (200 seconds "ON" followed by 400 seconds "OFF").

To assess the host specificity of the 23 phages, their infectivity was tested on four *Pseudomonas* spp. strains outside the *P. syringae* species complex using the double agar layer method on solid medium. The selected strains were *P. aeruginosa* CFBP 2466, *P. fluorescens* CFBP 2102, *P. graminis* 38B9, and *P. rhizosphaerae* 6B4. Strains were grown on King's B (KB) agar for 48 h at 25°C (*P. fluorescens*, *P. graminis*) or 28°C (*P. aeruginosa*, *P. rhizosphaerae*). Biological triplicates were inoculated into 5 mL of KB broth and incubated overnight at 25°C with shaking at 200 rpm. Then, 500 µL of each culture was transferred to 5 mL fresh KB and incubated under the same conditions for 4 h. 300 µL of the culture in exponential phase (OD between 0.2 and 0.8 at 600 nm) was mixed with 6 mL of 0.7% KB semi-solid agar for the overlay. Each phage at an initial concentration of $10^7$ PFU/mL was serially diluted, and 10 µL of each dilution were spotted in duplicate onto the bacterial lawn. Plates were allowed to dry and incubated overnight at 25°C before evaluation of plaque lysis. The phage concentrations on the tested strains were divided by the concentration obtained on their respective reference strains to calculate the efficiency of plating (EOP), a standard metric.

**b. Calculation of bacterial growth inhibition percentages and probability of significant inhibitory effects of phages on bacteria.** Bacterial growth curve data from each plate were analysed using R software (version 4.4.0) with the *growthcurver* package [87]. The area under the empirical growth curve (auce) was the primary parameter analysed. For each biological replicate of the bacterial control condition (without phage), the mean and standard deviation of the ace were calculated. The percentage of bacterial growth inhibition caused by each phage was then determined for each biological replicate using the following formula:
$$ \text{Bacterial growth inhibition}(\%) = \frac{\text{auce bacteria} - \text{auce phage}}{\text{auce bacteria}} \times 100 $$

where auce bacteria represent the mean auce of each bacterial control biological replicates, and auce phage corresponds to the auce of the same bacterial strain when exposed to a given phage. This normalisation accounts for potential differences in bacterial growth rates, allowing for direct comparison of phage effects across bacterial strains [88]. To evaluate potential bias in this calculation, particularly when bacteria exhibit different maximal densities [88], we assessed the coefficient of variation (CV) of the k parameter among bacteria control condition. The CV was low (9.22%), supporting the validity of this approach and indicating minimal bias associated with differences in growth maxima. Additionally, the CV of the auce values under bacteria control condition was moderate (11.34%), further supporting the robustness of the percentage inhibition metric for cross-strain comparisons while limiting potential biases [88].

Phage impact on bacterial growth was assessed using a probability value derived from a normalised distribution, calibrated based on the mean (μ in the formula) and standard deviation (σ in the formula) of the auce for each biological replicate of each bacterium exposed to phages. The probability ranges from 0 to 1 and is calculated using:

$$ \text{auce phage} \sim N(\mu(\text{auce bacteria}), \sigma(\text{auce bacteria})) $$

A probability value below 0.05 indicates a significant growth reduction in presence of phage, while a value exceeding 0.95 indicates a significant growth stimulation of bacteria in presence of phage. The distribution of inhibition percentages within the 0.05 to 0.95 probability range is approximately symmetric and unimodal (S8 Fig), resembling a Gaussian distribution. However, the Shapiro-Wilk test indicated a significant deviation from strict normality (p<0.01), likely due to the large sample size. Despite this, the distribution shape supports the use of this probability window to robustly summarise phage effects on bacterial growth. Bacterial strains 41A and 93D had only one technical replicate per biological replicate, precluding probability calculation.

A phage was considered capable of inhibiting a bacterial strain if at least one of the three bacterial replicates showed a significant inhibition (probability <0.05); in this case, the phage-bacterium combination was assigned a value of 1, in order to capture the maximal host range potential of each phage. If no replicates showed significant inhibition, the combination was assigned a value of 0. This threshold was verified to yield consistent results with a more stringent criterion (requiring at least two out of three replicates to be significant), in order to confirm the robustness of our analysis.

The relationship between the breadth of infection and the strength of inhibition was first examined using the full dataset. For the phage perspective, we plotted virulence (mean percentage of inhibition per phage) against infectivity (the number of bacterial strains each phage significantly infected) (S9A Fig). For the bacterial perspective, we plotted susceptibility (mean percentage of inhibition per bacterial strain) against sensitivity (the number of phages significantly infecting each strain) (S9B Fig). We then considered only combinations showing significant inhibition, for both phages (S9C Fig) and bacteria (S9D Fig), in order to reduce the dependency between variables. Spearman correlation tests were calculated exclusively on this subset of significant interactions, as this approach removes part of the inherent dependence generated by the shared underlying metric.

Generally, bacterial populations decreased in the presence of phages. In 17% of observed phage-bacteria combinations, a “stimulatory” effect was identified, with inhibition values reaching up to -40%. Given the lack of reproducibility of the stimulatory effect (S10 Fig), only positive inhibition values were retained for all analyses.

## Statistical analysis

**a. Analysis of the structure of phage-bacteria interaction networks.** The analysis of networks representing bacterial growth inhibition percentages of 44 bacterial strains combined with 23 phages was computed using the mean of
the three biological replicates for each combination. The analyses were performed using the R software (version 4.4.0) with the *bipartite* (for nestedness and modularity analysis) and *igraph* (for modularity analysis) packages, following Moury *et al.* [30] and Torres-Barceló *et al.* [46].

The nestedness of the inhibition percentage network was evaluated using two algorithms: *wNODF* (a weighted nestedness metric based on overlap and decreasing fill) [89] and *WINE* (weighted-interaction nestedness estimator) [90]. Modularity was assessed using six algorithms: *spinglass* [91–93], *edge betweenness* [91], *fast greedy* [94], *leading eigenvector* [95], *Beckett* [96], and *louvain* [97].

Prior to nestedness and modularity analyses, the continuous inhibitory percentage data of the matrix were transformed into discrete integer values spanning from 0 to 100.

In addition, we analysed a network of binary data (0 and 1 values) representing the statistical significance of bacterial growth inhibition by phages. For this network, three algorithms were used for nestedness: *WINE* as applied to quantitative data, *Binmatnest* [98], and *NODF2* [89] and three algorithms were used for modularity (*spinglass*, *leading eigenvector*, and *Beckett*).

For each algorithm and each network, the calculated nestedness and modularity scores were compared with the scores obtained with seven (network of continuous inhibition data) or six (network of binary inhibition significance data) different null models (i.e., simulated networks for which no particular nestedness or modularity patterns are expected) in order to estimate the statistical significance of nestedness and modularity of the actual networks [30,99,100]. Briefly, each null model differs by the constraints applied during randomization. The N model preserves the total sum and number of zeros but not row or column totals. C1 and R1 apply the same principle column- or row-wise, respectively. The B model maintains both row and column sums. The S model fully randomizes cell values, while C2 and R2 randomize values within columns or rows, respectively. The PD (Probabilistic Degree) model takes the original matrix and reshuffles positive infections while maintaining the same ‘degree distribution’, on average, of both phages and hosts while reshuffling. Together, these models span a gradient of constraint strength, from the least constrained (S, N) to the most constrained (B).

**b. Statistical analyses to explain bacterial sensitivity and susceptibility to phages.** All statistical analyses were conducted in R version 4.4.2. The phylogenetic signal of the number of phages infecting each bacterial strain was calculated using the *phylosig* function from the *phytools* package, based on a phylogenetic tree reconstructed from bacterial core-genome analysis [83].

The bacterial susceptibility profiles were tested with a PERmutational Multivariate ANalysis Of VAriance (PERMANOVA), a non-parametric method that assesses the influence of categorical variables on multivariate data by partitioning variation within a distance matrix. Unlike traditional ANOVA, it makes no assumptions about data distribution and relies on permutations to evaluate the statistical significance of group differences [101]. Here, Bray–Curtis dissimilarity was used to construct a distance matrix based on the mean percentage of bacterial growth inhibition (from three replicates per combination). This matrix was then analysed using the *adonis2* function from the *vegan* R package, with 9,999 permutations to ensure robust p-value estimation [102]. "Phylogroup", "ecological origin", and "whether the strain was an isolation bacterial strain for the phages" were included as explanatory variables. Multiple comparisons were conducted using the function *pairwise.adonis2* from the *pairwiseAdonis* package [103], which is specifically adapted for PERMANOVA post-hoc testing.

Additionally, the summed sensitivity, determined from the significant inhibition network (see above), was used as a response variable in a Generalised Linear Model (GLM) with a Poisson distribution with the function *glmmTMB* from the package *glmmTMB* [104]. The model included "phylogroup", "ecological origin", and their interaction as fixed effects. As no significant interaction was detected, the simplified additive model (*phylogroup* + *ecological origin*) was retained. No overdispersion was observed in the data (deviance/degrees of freedom = 1.01). The effects of each explanatory variable were assessed using Type III ANOVA with the *Anova* function from the *car* package. Post-hoc pairwise comparisons were performed using Tukey's test implemented via the *multcomp* and *multcompView* packages, with significance groups assigned accordingly.
**c. Statistical analysis to explain infectivity and virulence of phages.** The Phylogenetic Host Range Index (PHRI) [46], derived from the phylogenetic diversity index *Hp* [105], was calculated using the *hilldiv* package. For this analysis, the bacterial phylogenetic tree, constructed from core-genome analysis (see above), was transformed into an ultrametric tree, and the PHRI index was determined using the matrix of percentage inhibition.

The ranking of phage genera based on inhibition percentage (offset by +0.1 to avoid zeros) was assessed using a GLMM with a Gamma distribution with a log-link, suitable for strictly positive, continuous, and skewed data. The model included a random effect accounting for the "phage" tested, the "bacterial" strain tested, and "replicates" nested within each "phage:bacterium" combination. No overdispersion was detected (χ²/ddl=0.67). Post-hoc pairwise comparisons between genera were performed using Tukey's test.

The phage-host infection profiles were tested with a PERMANOVA test based on the distance matrix based on the mean percentage of bacterial growth inhibition compared with "phage genus" and "whether the phage had been isolated using one of the bacterial strains included in the experimental panel" explanatory variables (Bray-Curtis distance matrix; 9,999 permutations). Pairwise comparisons of phage infection patterns were performed using the *pairwise.adonis2* function [103].

**d. Statistical analysis of defence systems and prophages.** The phylogenetic signal of the total number of defence systems and the abundance of prophages was calculated using the *phylosig* function from the *phytools* package, based on a phylogenetic tree reconstructed from bacterial core-genome analysis [83].

The total number of bacterial defence systems per bacterial genome and the abundance of prophages were tested against the explanatory variables "phylogroup", "ecological origin", and "growth rate" in a GLM with a Poisson distribution with an additive model of the three explanatory variables. No overdispersion was observed in the data (deviance/degrees of freedom for defence systems = 0.61; for abundance of prophages = 0.47), and the effects of each explanatory variable were assessed using Type III ANOVA. Finally, the post-hoc pairwise comparisons were performed using Tukey's test and significance groups were assigned accordingly.

The relationship between ANI of bacterial pairs and their percentage similarity in their defence systems (offset by +0.1 to avoid zeros) was assessed using a GLMM with a Gamma distribution and a log-link function. Random effects for each bacterium in the pair were included to account for the non-independence of observations. No overdispersion was detected (χ²/ddl=0.24).

The defence system profiles, the proportions of the defence system categories, and prophages distribution profile across bacterial strains were analysed using a PERMANOVA test based on a Bray-Curtis dissimilarity matrix with 9,999 permutations. "Phylogroup", "ecological origin", and "growth rate" were included as explanatory variables. Pairwise comparisons of the defence system and prophage composition between groups were conducted using the *pairwise.adonis2* function [103].

Linear regressions were conducted to explain the number of phages infecting each bacterial strain, using the *total number of defence system subtypes per bacterium*, the *phylogroups*, and *ecological origins* as explanatory variables in a multiple linear regression model (adjusted R²=0.21; p=0.02 for the overall model). The same analysis was applied using the *number of prophages per bacterium* instead of defence systems, with *phylogroups* and *ecological origins* as predictors (adjusted R²=-0.06; p=0.81 for the overall model). Before conducting these analyses, all assumptions of linear regression, including normality of residuals and homoscedasticity, were verified.

**e. Statistical analysis of host range data with multimodel analysis.** To account for the widest range of explanatory variables influencing the percentage of inhibition, we performed a multimodel analysis coupled with a Generalised Linear Mixed Model (GLMM). To enhance model efficiency, the number of explanatory variables was reduced through a multi-step process. Variables with more than 20% missing data were excluded. To ensure independence among explanatory variables before multimodel analysis, we assessed pairwise associations using Spearman's correlation for quantitative variables, chi-squared tests for categorical variables, and Kruskal-Wallis tests for associations between categorical and quantitative
<figure><img src="image_17.png" /></figure>

variables. Variables showing significant correlations (p<0.05) were considered dependent. In such cases, one of the correlated variables was retained, while the others were excluded. To further refine the selection, Principal Component Analyses (PCA) were conducted separately on the matrices representing the distributions of defence systems and prophages across bacterial genomes. The most informative variables contributing to the first principal components were retained for subsequent modelling.

The PCA on defence systems was performed without centring or scaling the data. Two dimensions accounted for 32.5% of the variance and led to the selection of the variables "Septu", "Gabija", and "RloC", which were among the most frequent defence systems in the genomes. Similarly, the PCA on the distribution of prophage species in bacterial genomes, also conducted without centring or scaling (already binary), captured 14% of the variance in the first two dimensions and identified two perpendicular variable groups composed respectively of prophage species "17", "54", "26", and "29" on one side, and species "7", "47", "57", "22", and "31" on the other. One representative species from each group was selected for the multimodel analysis (species "17" and "22").

To identify factors explaining variations in bacterial growth inhibition percentage, we fitted a GLMM. The response variable (offset by +0.1 to avoid zeros) was modelled using a Gamma distribution with a log-link function. Random effects included "phage", "bacterium", and their "imbricated replicates".

To compare two approaches for modelling the effects of bacterial defence systems and prophages on the percentage of inhibition, two separate GLMMs were constructed: one using total counts of defence systems and prophages, and another incorporating PCA-derived principal components of the distribution matrices (of prophages and defence systems). No overdispersion was detected in either model (total counts: $\chi^2/df = 0.65$; PCA-derived: $\chi^2/df = 0.66$), supporting the appropriateness of the Gamma distribution for modelling the response variable.

To identify the most informative predictors explaining variation in phage replication efficiency, we used a multimodel analysis. Model selection was performed using the *dredge* function from the *MuMln* package, ranking models based on corrected Akaike information criterion (AIC). All possible sub-models were generated and evaluated. Models whose cumulative AIC weights accounted for 95% of the total were retained for model averaging using the *model.avg* function. The relative importance of each explanatory variable was then estimated using the *sw* function, which computes the sum of the AIC weights across all models in which the variable appears.

## Supporting information

**S1 Fig. Bacterial resistance to phages and host range breadth of the tested phages.** (A) Number of phages to which bacteria are sensitive or resistant, with the phylogeny of the strains indicated at the bottom of the histogram. (B) Number of bacteria infected or not by the phages, with the taxonomy of the phages listed at the bottom of the histogram. Values are derived from the binary inhibition matrix.

(JPG)

**S2 Fig. Taxonomic and phylogenetic characteristics, and ecological origins of phages and bacteria explaining host range.** (A) Bacterial percentage of inhibition by phages across bacterial phylogroups and ecological origins. Significance letters indicate results from the PERMANOVA test and associated pairwise comparisons. (B) Bacterial inhibition percentage of different phage genera. Significance letters represent results from a GLMM with a Gamma regression distribution and random effect (phage, bacteria, and the replicate imbricated in each of the phage-bacteria interactions), and post hoc pairwise comparisons performed using Tukey's test.

(JPG)

**S3 Fig. Relative proportions of defence systems by functional category (%) and total number of defence systems across bacterial genomes.** On the left, a phylogenetic tree based on the core genomes of the bacterial strains is shown, with corresponding phylogroups indicated in the adjacent column.

(JPG)
**S4 Fig.** Percentage of similarity of defence system subtypes among bacterial strains in relation to their genome-wide average nucleotide identity (ANI). Transparent grey dots represent individual observations (darker areas indicate overlap). The overall trend is shown by the grey line, which represents the mean prediction from the GLMM accounting for random effects associated with bacterial strains. The shaded ribbon indicates the 95% confidence interval of the model's marginal predictions.

(JPG)

**S5 Fig.** Co-phylogeny of the CL4 prophage family-like (R-type syringacin). Co-phylogenetic analysis between 68 strains of the *P. syringae* species complex (left), represented by a maximum likelihood phylogenetic tree based on their core-genomes, and their associated 68 CL4 family-like prophages (right), represented by a proteomic tree. The 44 bacterial strains included in this study are marked with dots at the tips of the bacterial tree. Phylogroups and prophage-carrying strains are colour-coded accordingly. Two nodes with a maximum likelihood support value below 0.90 are marked with red dots.

(JPG)

**S6 Fig.** Experimental setup illustrating bacterial growth inhibition by phages. Each biological replicate per bacterial strain is represented by a different colour, and each row corresponds to the technical replicates within that biological replicate. For clarity, only 3 of the 23 tested phages are shown in this schematic, which does not depict the randomisation of phages across plates. Three types of internal controls are included: positive controls with bacteria only, negative controls with medium only, and internal reference controls using the bacterial strain P1.01.01B03 (growth shown in S7B Fig) and the internal phage Ghuch01, which consistently inhibits this strain (see S7C Fig). Icons were obtained from open-source resources (microplate: https://commons.wikimedia.org/wiki/File:96-Well Plate.svg; phage: https://openclipart.org/detail/62785/virus; bacteria: https://openclipart.org/detail/221334/lamarque-disease; bottle of media: https://openclipart.org/detail/128449/reagent-bottle).

(JPG)

**S7 Fig.** Analysis of the different controls used in the phage host range *in vitro* tests. (A) Analysis of culture medium controls by plate (labelled by bacteria). (B) Area under the empirical curve (auc) of the internal bacterial control, *P. syringae* strain P1.01.01B03. (C) Inhibition of the internal phage control Ghuch01 on the internal bacterial control (%).

(JPG)

**S8 Fig.** Histogram showing the distribution of inhibition percentages within the 0.05 to 0.95 probability range, displaying an approximately symmetric and unimodal shape.

(JPG)

**S9 Fig.** Relationship between infection breadth of infection and inhibition strength. (A, C) Phage perspective: mean percentage of inhibition per phage (y-axis) is plotted against the number of bacterial strains each phage significantly infected (x-axis). Colours indicate phage genera. (B, D) Bacterial perspective: mean percentage of inhibition per bacterial strain (y-axis) is plotted against the number of phages significantly infecting each strain (x-axis). Colours represent bacterial phylogroups, and shapes indicate ecological niches. (A, B) Analyses based on overall inhibition values. (C, D) Analyses restricted to significant phage-bacterium interactions. Spearman's rho and associated p-values are reported for panels (C) and (D). Icons were obtained from open-source resources (phage: https://openclipart.org/detail/62785/virus; bacteria: https://openclipart.org/detail/221334/lamarque-disease).

(JPG)

**S10 Fig.** Results of additional experiments on five potentially stimulating phages and five potentially stimulated bacteria. Generally, bacterial populations decreased in the presence of phages. In 17% of observed
phage-bacteria combinations, a “stimulatory” effect was identified, with inhibition values reaching up to -40%. Five phages showing strong stimulatory effect and five bacteria showing significant stimulation were chosen to test the reproducibility of the stimulation effect in subsequent experiments. (A) Matrix of mean inhibition percentages for the five *P. syringae* strains (rows) and five phages (columns), with their phylogenetic and taxonomic affiliations indicated by colour. The number of stars indicates the number of replicates showing a significantly inhibiting or stimulating percentage. Of ten initially identified significant stimulating interactions (delimited by red lines), only three were confirmed in subsequent experiments. (B) Size of the five bacterial populations in response to the phages. Stars correspond to significant differences from the control (no phage) based on the results of the GLM model, family = Gamma. Arrows indicate interactions with at least one of the replicates showing a significant inhibition or stimulation percentage. Out of the three significant interactions of stimulatory phenotype, only 2 phage-bacteria interactions resulted in slightly higher bacterial densities in the presence of phages compared to phage-free conditions, though this difference was not statistically significant. (C) Spearman correlation between bacterial concentration and inhibition percentage, by bacterium and phage. No correlation was significant and bacterial population size only poorly explains the inhibition percentages.

(JPG)

**S11 Fig. Kinetics of bacterial growth inhibition across the full dataset.** Optical density (OD) is plotted over time (hours) for each phage-bacteria interaction, with phages (n=23) in columns and bacterial strains (n=44) in rows. Grey curves represent the average growth of bacteria alone for each phage-bacteria pair. Blue curves show individual replicates for each interaction (n=3), illustrating the variability in bacterial response to phage infection.

(JPG)

**S1 Table. Phage-bacterium matrices showing (A) the percentage of bacterial growth inhibition and (B) the number of biological replicates (out of three) showing significant inhibition.** Rows correspond to bacterial strains, while columns represent phage clones, each annotated with their respective phylogroups and genera.

(XLSX)

**S2 Table. Analysis of nestedness and modularity of the network of binary data representing the significance or non-significance of bacterial growth inhibition by phages (S1B Table).** Null models are described in Moury *et al.* [30] for B, N, C2, R2 and S, and in Moury *et al.* [99] for PD. For each null model and each algorithm, the frequency of simulations, among 100, showing a strictly higher nestedness (or modularity) score than that of the actual matrix is indicated. Significantly nested or modular networks ($p \le 0.05$) are indicated in red.

(XLSX)

**S3 Table. Genetic analysis of the 88 prophages identified in the 44 *P. syringae* strains using VIRIDIC.** The 27 reference phages are highlighted in yellow (25 from the collection described in Feltin *et al.* [72] and 3 from NCBI). Species and genus-level clusters of the prophages are indicated on the right in colour. High genetic similarities are highlighted in green.

(XLSX)

**S4 Table. Phylogenetic characteristics, ecological origins, countries, and dates of sampling of the 44 *P. syringae* strains used for the host range study.** The four *P. syringae* strains highlighted in bold are strains that were used to isolate some of the phages tested in the host range assay. Four strains are in the CIRM-CFBP, the French collection for plant-associated bacteria (https://cirm-cfbp.fr/; June 2025): USA0007(CFBP 8574), CC1559 (CFBP 8551), MAFF302273PT (CFBP 2339), and CC1583 (CFBP 8564).

(XLSX)
**Acknowledgments**

We thank Elisabeth Villard, Cécile Monteil, and Charlotte Chandeysson for their assistance with experimental procedures.
We are also grateful to Dominique Holtappels for valuable discussions, and to Hannah Toth for her assistance in deposit-
ing bacterial genomes into NCBI.

**Author contributions**

**Conceptualization:** Chloé Feltin, Cindy E. Morris, Clara Torres-Barceló.

**Data curation:** Chloé Feltin, Sylvain Piry, Lola Chateau, Theo H. M. Smits.

**Formal analysis:** Chloé Feltin, Sylvain Piry, Benoit Moury, Lola Chateau, Karine Berthier, Jonathan M. Jacobs, Jules Butchacas, Linda Fenske, Lillian Ebeling-Koning, Theo H. M. Smits, Clara Torres-Barceló.

**Funding acquisition:** Cindy E. Morris, Clara Torres-Barceló.

**Investigation:** Chloé Feltin, Sylvain Piry, Benoit Moury, Lola Chateau, Karine Berthier, Jules Butchacas, Linda Fenske, Lillian Ebeling-Koning, Theo H. M. Smits, Cindy E. Morris, Clara Torres-Barceló.

**Methodology:** Chloé Feltin, Lola Chateau, Cindy E. Morris, Clara Torres-Barceló.

**Project administration:** Chloé Feltin, Clara Torres-Barceló.

**Resources:** Clara Torres-Barceló.

**Supervision:** Clara Torres-Barceló.

**Writing – original draft:** Chloé Feltin, Clara Torres-Barceló.

**Writing – review & editing:** Chloé Feltin, Sylvain Piry, Benoit Moury, Jonathan M. Jacobs, Theo H. M. Smits, Cindy E. Morris, Clara Torres-Barceló.

**References**

1. Brown TL, Charity OJ, Adriaenssens EM. Ecological and functional roles of bacteriophages in contrasting environments: marine, terrestrial and human gut. Curr Opin Microbiol. 2022;70:102229.

2. Wilhelm SW, Suttle CA. Viruses and nutrient cycles in the sea: viruses play critical roles in the structure and function of aquatic food webs. BioScience. 1999;49(10):781–8.

3. Wommack KE, Colwell RR. Virioplankton: viruses in aquatic ecosystems. Microbiol Mol Biol Rev. 2000;64(1):69–114. https://doi.org/10.1128/MMBR.64.1.69-114.2000 PMID: 10704475

4. Suttle CA. Marine viruses--major players in the global ecosystem. Nat Rev Microbiol. 2007;5(10):801–12. https://doi.org/10.1038/nrmicro1750 PMID: 17853907

5. Brussaard CPD, Wilhelm SW, Thingstad F, Weinbauer MG, Bratbak G, Heldal M, et al. Global-scale processes with a nanoscale drive: the role of marine viruses. ISME J. 2008;2(6):575–8. https://doi.org/10.1038/ismej.2008.31 PMID: 18385772

6. Mirzaei MK, Maurice CF. Ménage à trois in the human gut: interactions between host, bacteria and phages. Nat Rev Microbiol. 2017;15(7):397–408. https://doi.org/10.1038/nrmicro.2017.30 PMID: 28461690

7. De Sordi L, Lourenço M, Debarbieux L. The Battle Within: Interactions of Bacteriophages and Bacteria in the Gastrointestinal Tract. Cell Host Microbe. 2019;25(2):210–8. https://doi.org/10.1016/j.chom.2019.01.018 PMID: 30763535

8. Svircev A, Roach D, Castle A. Framing the Future with Bacteriophages in Agriculture. Viruses. 2018;10(5):218. https://doi.org/10.3390/v10050218 PMID: 29693561

9. Hyman P, Abedon ST. Bacteriophage host range and bacterial resistance. In: Advances in Applied Microbiology [Internet]. Elsevier; 2010 [cited 2024 Aug 12]. p. 217–48. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0065216410700071

10. Georjon H, Bernheim A. The highly diverse antiphage defence systems of bacteria. Nat Rev Microbiol. 2023;21(10):686–700. https://doi.org/10.1038/s41579-023-00934-x PMID: 37460672

11. Piel D, Bruto M, Labreuche Y, Blanquart F, Goudenège D, Barcia-Cruz R, et al. Phage-host coevolution in natural populations. Nat Microbiol. 2022;7(7):1075–86. https://doi.org/10.1038/s41564-022-01157-1 PMID: 35760840
12. Borin JM, Lee JJ, Lucia-Sanz A, Gerbino KR, Weitz JS, Meyer JR. Rapid bacteria-phage coevolution drives the emergence of multiscale networks. Science. 2023;382(6671):674-8.

13. Beamud B, García-González N, Gómez-Ortega M, González-Candelas F, Domingo-Calap P, Sanjuan R. Genetic determinants of host tropism in Klebsiella phages. Cell Rep. 2023;42(2):112048. https://doi.org/10.1016/j.celrep.2023.112048 PMID: 36753420

14. Gaborieau B, Vaysset H, Tesson F, Charachon I, Dib N, Bernier J, et al. Prediction of strain level phage-host interactions across the Escherichia genus using only genomic information. Nat Microbiol. 2024;9(11):2847-61. https://doi.org/10.1038/s41564-024-01832-5 PMID: 39482383

15. Johnson MC, Laderman E, Huiting E, Zhang C, Davidson A, Bondy-Denomy J. Core defense hotspots within Pseudomonas aeruginosa are a consistent and rich source of anti-phage defense systems. Nucleic Acids Res. 2023;51(10):4995-5005. https://doi.org/10.1093/nar/gkad317 PMID: 37140042

16. Burke KA, Urick CD, Mzhavia N, Nikolich MP, Filippov AA. Correlation of Pseudomonas aeruginosa Phage Resistance with the Numbers and Types of Antiphage Systems. Int J Mol Sci. 2024;25(3):1424. https://doi.org/10.3390/ijms25031424 PMID: 38338703

17. Costa AR, van den Berg DF, Esser JQ, Muralidharan A, van den Bossche H, Bonilla BE, et al. Accumulation of defense systems in phage-resistant strains of Pseudomonas aeruginosa. Sci Adv. 2024;10(8):eadj0341. https://doi.org/10.1126/sciadv adj0341 PMID: 38394193

18. Hussain FA, Dubert J, Elsherbini J, Murphy M, VanInsberghe D, Arevalo P. Rapid evolutionary turnover of mobile genetic elements drives bacterial resistance to phages. Science. 2021;374(6566):488-92.

19. Czajkowski R. May the phage be with you? Prophage-like elements in the genomes of soft rot Pectobacteriaceae: Pectobacterium spp. and Dickeya spp. Front Microbiol. 2019. https://doi.org/10.3389/fmicb.2019.00138

20. Dominguez-Mirazo M, Jin R, Weitz JS. Functional and Comparative Genomic Analysis of Integrated Prophage-Like Sequences in "Candidatus Liberibacter asiaticus". mSphere. 2019;4(6):e00409-19. https://doi.org/10.1128/mSphere.00409-19 PMID: 31722990

21. Castillo JA, Secaira-Morocho H, Maldonado S, Sarmiento KN. Diversity and evolutionary dynamics of antiphage defense systems in *Ralstonia solanacearum* species complex. Front Microbiol [Internet]. 2020 [cited 2024 Sept 24];11. Available from: https://www.frontiersin.org/journals/ microbiology/articles/10.3389/fmicb.2020.00961/full

22. Greenrod STE, Stoycheva M, Elphinstone J, Friman V-P. Global diversity and distribution of prophages are lineage-specific within the Ralstonia solanacearum species complex. BMC Genomics. 2022;23(1):689. https://doi.org/10.1186/s12864-022-08909-7 PMID: 36199029

23. Dougherty PE, Nielsen TK, Riber L, Lading HH, Forero-Junco LM, Kot W, et al. Widespread and largely unknown prophage activity, diversity, and function in two genera of wheat phyllosphere bacteria. ISME J. 2023;17(12):2415-25. https://doi.org/10.1038/s41396-023-01547-1 PMID: 37919394

24. Holtappels D, Rickus GEJ, Morgan T, de Rezende RR, Koskella B, Alfenas-Zerbini P. Evidence for the key roles of the *Pseudomonas syringae* mobilome in shaping biotic interactions [Internet]. bioRxiv; 2024 [cited 2024 Aug 20]. p. 2024.03.19.585818. Available from: https://www.biorxiv.org/content/10.1101/2024.03.19.585818v2

25. Jimenez D, Beltran D, Castillo JA. Mobile genetic elements of Xylella fastidiosa and their contribution to pathogenicity. Plant Pathol. 2024;73(9):2490-9. https://doi.org/10.1111/ppa.13992

26. Flores CO, Meyer JR, Valverde S, Farr L, Weitz JS. Statistical structure of host-phage interactions. Proc Natl Acad Sci U S A. 2011;108(28):E288-97. https://doi.org/10.1073/pnas.1101595108 PMID: 21709225

27. Gupta A, Peng S, Leung CY, Borin JM, Medina SJ, Weitz JS, et al. Leapfrog dynamics in phage-bacteria coevolution revealed by joint analysis of cross-infection phenotypes and whole genome sequencing. Ecol Lett. 2022;25(4):876-88. https://doi.org/10.1111/ele.13965 PMID: 35092147

28. Molina F, Menor-Flores M, Fernández L, Vega-Rodríguez MA, García P. Systematic analysis of putative phage-phage interactions on minimum-sized phage cocktails. Sci Rep. 2022;12(1):2458. https://doi.org/10.1038/s41598-022-06422-1 PMID: 35165352

29. Kauffman KM, Chang WK, Brown JM, Hussain FA, Yang J, Polz MF, et al. Resolving the structure of phage-bacteria interactions in the context of natural diversity. Nat Commun. 2022;13(1):372. https://doi.org/10.1038/s41467-021-27583-z PMID: 35042853

30. Moury B, Audergon JM, Baudracco-Arnas S, Ben Krima S, Bertrand F, Boissot N. The quasi-universality of nestedness in the structure of quantitative plant-parasite interactions. Peer Community J. 2021;1:e44.

31. Morris CE, Monteil CL, Berge O. The life history of Pseudomonas syringae: linking agriculture to earth system processes. Annu Rev Phytopathol. 2013;51:85-104. https://doi.org/10.1146/annurev-phyto-082712-102402 PMID: 23663005

32. Morris CE, Sands DC, Vinatzer BA, Glaux C, Guilbaud C, Buffière A, et al. The life history of the plant pathogen Pseudomonas syringae is linked to the water cycle. ISME J. 2008;2(3):321-34. https://doi.org/10.1038/ismej.2007.113 PMID: 18185595

33. Morris CE, Lamichhane JR, Nikolić I, Stanković S, Moury B. The overlapping continuum of host range among strains in the Pseudomonas syringae complex. Phytopathol Res. 2019;1(1):4.

34. Parisi L, Morgaint B, Blanco-Garcia J, Guilbaud C, Chandeysson C, Bourgeay JF, et al. Bacteria from four phylogroups of the Pseudomonas syringae complex can cause bacterial canker of apricot. Plant Pathol. 2019;68(7):1249-58. https://doi.org/10.1111/ppa.13051

35. Thomas MD. Bacteriophages from Sewage Specific for Fluorescent Phytopathogenic Pseudomonads. Phytopathology. 1983;73(3):403. https://doi.org/10.1094/phyto-73-403

36. Cuppels DA. Isolation and Characterization of Phages Useful for Identifying Pseudomonas syringaepv.tomato. Phytopathology. 1983;73(10):1376. https://doi.org/10.1094/phyto-73-1376
37. Gross DC, Cody YS, Proebsting EL, Radamaker GK, Spotts RA. Ecotypes and pathogenicity of ice-nucleation-active Pseudomonas syringae isolated from deciduous fruit tree orchards. Phytopathology. 1983;74:241-8.

38. Cuppels DA. The use of pathovar-indicative bacteriophages for rapidly detecting Pseudomonas syringae pv. tomato in tomato leaf and fruit lesions. Phytopathology. 1984;74(2):891-4.

39. Taylor JD, Teverson DM, Allen DJ, Pastor-Corrales MA. Identification and origin of races of Pseudomonas syringae pv. phaseolicola from Africa and other bean growing areas. Plant Pathol. 1996;45(3):469-78.

40. Frampton RA, Taylor C, Holguín Moreno AV, Visnovsky SB, Petty NK, Pitman AR, et al. Identification of bacteriophages for biocontrol of the kiwifruit canker phytopathogen Pseudomonas syringae pv. actinidiae. Appl Environ Microbiol. 2014;80(7):2216-28. https://doi.org/10.1128/AEM.00062-14 PMID: 24487530

41. James SL, Rabiey M, Neuman BW, Percival G, Jackson RW. Isolation, characterisation and experimental evolution of phage that infect the horse chestnut tree pathogen, Pseudomonas syringae pv. aesculi. Curr Microbiol. 2020;77(8):1438-47.

42. Rombouts S, Volckaert A, Venneman S, Declercq B, Vandenhevel D, Allonsius CN, et al. Characterization of novel bacteriophages for biocontrol of bacterial blight in leek caused by *Pseudomonas syringae* pv. *porri*. Front Microbiol. 2016.

43. Geng Y, Nguyen TVP, Homaee E, Golding I. Using bacterial population dynamics to count phages and their lysogens. Nat Commun. 2024;15(1):7814. https://doi.org/10.1038/s41467-024-51913-6 PMID: 39242585

44. Wright RCT, Friman V-P, Smith MCM, Brockhurst MA. Cross-resistance is modular in bacteria-phage interactions. PLoS Biol. 2018;16(10):e2006057. https://doi.org/10.1371/journal.pbio.2006057 PMID: 30281587

45. Wright RCT, Friman V-P, Smith MCM, Brockhurst MA. Resistance Evolution against Phage Combinations Depends on the Timing and Order of Exposure. mBio. 2019;10(5):e01652-19. https://doi.org/10.1128/mBio.01652-19 PMID: 31551330

46. Torres-Barceló C, Boyer C, Garneau JR, Poussier S, Robène I, Moury B. A phylogenetic host-range index reveals ecological constraints in phage specialisation and virulence. Mol Ecol. 2025;34(18):e70052. https://doi.org/10.1111/mec.70052

47. Seidel R, Bloom JGP, Dekker C, Szczelkun MD. Motor step size and ATP coupling efficiency of the dsDNA translocase EcoR124I. EMBO J. 2008;27(9):1388-98. https://doi.org/10.1038/emboj.2008.69 PMID: 18388857

48. Pleška M, Qian L, Okura R, Bergmiller T, Wakamoto Y, Kussell E. Bacterial autoimmunity due to a restriction-modification system. Curr Biol. 2016;26(3):404-9.

49. Berngruber TW, Lion S, Gandon S. Evolution of suicide as a defence strategy against pathogens in a spatially structured environment. Ecol Lett. 2013;16(4):446-53. https://doi.org/10.1111/ele.12064 PMID: 23331662

50. van Houte S, Buckling A, Westra ER. Evolutionary ecology of prokaryotic immune mechanisms. Microbiol Mol Biol Rev. 2016;80(3):745-63.

51. Turner D, Kropinski AM, Adriaenssens EM. A Roadmap for Genome-Based Phage Taxonomy. Viruses. 2021;13(3):506. https://doi.org/10.3390/v13030506 PMID: 33803862

52. Mavrodi DV, Loper JE, Paulsen IT, Thomashow LS. Mobile genetic elements in the genome of the beneficial rhizobacterium Pseudomonas fluorescens Pf-5. BMC Microbiol. 2009;9(1):8.

53. Hockett KL, Renner T, Baltrus DA. Independent Co-Option of a Tailed Bacteriophage into a Killing Complex in Pseudomonas. mBio. 2015;6(4):e00452. https://doi.org/10.1128/mBio.00452-15 PMID: 26265717

54. Backman T, Latorre SM, Symeonidi E, Muszyński A, Bleak E, Eads L, et al. A phage tail-like bacteriocin suppresses competitors in metapopulations of pathogenic bacteria. Science. 2024;384(6701):eado0713. https://doi.org/10.1126/science.ado0713 PMID: 38870284

55. Fautt C, Hockett K, Delattre S, Baltrus D, Couradeau E. Tailocin tail fiber diversity correlates with rfbD variation in the Pseudomonas syringae species complex. ISME Commun. 2025;5(1):ycaf099. https://doi.org/10.1093/ismeco/ycaf099 PMID: 40589794

56. Hunter M, Fusco D. Superinfection exclusion: A viral strategy with short-term benefits and long-term drawbacks. PLoS Comput Biol. 2022;18(5):e1010125. https://doi.org/10.1371/journal.pcbi.1010125 PMID: 35536864

57. Koskella B, Lin DM, Buckling A, Thompson JN. The costs of evolving resistance in heterogeneous parasite environments. Proc R Soc B Biol Sci. 2011;279(1735):1896-903.

58. Morris CE, Sands DC, Vanneste JL, Montarry J, Oakley B, Guilbaud C. Inferring the evolutionary history of the plant pathogen Pseudomonas syringae from its biogeography in headwaters of rivers in North America, Europe, and New Zealand. mBio. 2010;1(3):e00107-10. https://doi.org/10.1128/mBio.00107-10

59. Monteil CL, Yahara K, Studholme DJ, Mageiros L, Méric G, Swingle B, et al. Population-genomic insights into emergence, crop adaptation and dissemination of Pseudomonas syringae pathogens. Microb Genom. 2016;2(10):e000089. https://doi.org/10.1099/mgen.0.000089 PMID: 28348830

60. Berge O, Monteil CL, Bartoli C, Chandeysson C, Guilbaud C, Sands DC, et al. A user's guide to a data base of the diversity of Pseudomonas syringae and its application to classifying strains in this phylogenetic complex. PLoS One. 2014;9(9):e105547. https://doi.org/10.1371/journal.pone.0105547 PMID: 25184292

61. Koonin EV, Makarova KS, Wolf YI. Evolutionary genomics of defense systems in archaea and bacteria. Annu Rev Microbiol. 2017;71(71):233-61.

62. Bernheim A, Sorek R. The pan-immune system of bacteria: antiviral defence as a community resource. Nat Rev Microbiol. 2020;18(2):113-9. https://doi.org/10.1038/s41579-019-0278-2 PMID: 31695182
63. Makarova KS, Wolf YI, Koonin EV. Comparative genomics of defense systems in archaea and bacteria. Nucleic Acids Res. 2013;41(8):4360-77.
https://doi.org/10.1093/nar/gkt157 PMID: 23470997

64. Cooney-Nutley K, Chakravorty S, Nix I, Zeng Z, Greer SF, Rabiey M. Understanding Bacterial Antiviral Defence Systems and Phage Receptors to Better Inform Rational Phage Cocktail Design to Treat Bacterial Canker. Microb Biotechnol. 2025;18(9):e70232. https://doi.org/10.1111/1751-7915.70232 PMID: 40970513

65. Xie Y, Wahab L, Gill JJ. Development and Validation of a Microtiter Plate-Based Assay for Determination of Bacteriophage Host Range and Virulence. Viruses. 2018;10(4):189. https://doi.org/10.3390/v10040189 PMID: 29649135

66. Bell G. Evolutionary rescue. Annu Rev Evol Syst. 2017;48(2017):605-27.

67. Abedon ST. Lysis from without. Bacteriophage. 2011;1(1):46-9. https://doi.org/10.4161/bact.1.1.13980 PMID: 21687534

68. Nobrega FL, Vlot M, de Jonge PA, Dreesens LL, Beaumont HJE, Lavigne R, et al. Targeting mechanisms of tailed bacteriophages. Nat Rev Microbiol. 2018;16(12):760-73. https://doi.org/10.1038/s41579-018-0070-8 PMID: 30104690

69. Bertozzi Silva J, Storms Z, Sauvageau D. Host receptors for bacteriophage adsorption. FEMS Microbiol Lett. 2016;363(4):fnw002. https://doi.org/10.1093/femsle/fnw002 PMID: 26755501

70. Müller DM, Pourtois JD, Kim MK, Targ B, Burgener EB, Milla C, et al. Bacterial receptors but not anti-phage defence mechanisms determine host range for a pair of *Pseudomonas aeruginosa* lytic phages [Internet]. bioRxiv; 2024 [cited 2024 Sept 23]. p. 2024.04.30.591980. Available from:
https://www.biorxiv.org/content/10.1101/2024.04.30.591980v2

71. Kering KK, Kibii BJ, Wei H. Biocontrol of phytobacteria with bacteriophage cocktails. Pest Manag Sci. 2019;75(7):1775-81. https://doi.org/10.1002/ps.5324 PMID: 30624034

72. Feltin C, Garneau JR, Morris CE, Bérard A, Torres-Barceló C. Novel phages of Pseudomonas syringae unveil numerous potential auxiliary metabolic genes. J Gen Virol. 2024;105(6):001990. https://doi.org/10.1099/jgv.0.001990 PMID: 38833289

73. King EO, Ward MK, Raney DE. Two simple media for the demonstration of pyocyanin and fluorescin. J Lab Clin Med. 1954;44(2):301-7.

74. Dieckmann MA, Beyvers S, Nkouamedjo-Fankep RC, Hanel PHG, Jelonek L, Blom J. EDGAR3.0: comparative genomics and phylogenomics on a scalable infrastructure. Nucleic Acids Res. 2021;49(W1):W185-92.

75. Price MN, Dehal PS, Arkin AP. FastTree 2--approximately maximum-likelihood trees for large alignments. PLoS One. 2010;5(3):e9490. https://doi.org/10.1371/journal.pone.0009490 PMID: 20224823

76. Tesson F, Hervé A, Mordret E, Touchon M, d'Humières C, Cury J, et al. Systematic and quantitative view of the antiviral arsenal of prokaryotes. Nat Commun. 2022;13(1):2561. https://doi.org/10.1038/s41467-022-30269-9 PMID: 35538097

77. Camargo AP, Roux S, Schulz F, Babinski M, Xu Y, Hu B. Identification of mobile genetic elements with geNomad. Nat Biotechnol. 2023;1-10.

78. Nayfach S, Camargo AP, Schulz F, Eloe-Fadrosh E, Roux S, Kyrpides NC. CheckV assesses the quality and completeness of metagenome-assembled viral genomes. Nat Biotechnol. 2021;39(5):578-85. https://doi.org/10.1038/s41587-020-00774-7 PMID: 33349699

79. Moraru C, Varsani A, Kropinski AM. VIRIDIC-A Novel Tool to Calculate the Intergenomic Similarities of Prokaryote-Infecting Viruses. Viruses. 2020;12(11):1268. https://doi.org/10.3390/v12111268 PMID: 33172115

80. Nishimura Y, Yoshida T, Kuronishi M, Uehara H, Ogata H, Goto S. ViPTree: the viral proteomic tree server. Bioinforma Oxf Engl. 2017;33(15):2379-80.

81. Orfei B, Pothier JF, Fenske L, Blom J, Moretti C, Buonaurio R, et al. Race-specific genotypes of Pseudomonas syringae pv. tomato are defined by the presence of mobile DNA elements within the genome. Front Plant Sci. 2023;14:1197706. https://doi.org/10.3389/fpls.2023.1197706 PMID: 37476164

82. Balbuena JA, Míguez-Lozano R, Blasco-Costa I. PACo: a novel procrustes application to cophylogenetic analysis. PLoS One. 2013;8(4):e61048.
https://doi.org/10.1371/journal.pone.0061048 PMID: 23580325

83. Revell LJ. Phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and other things). PeerJ. 2024;12:e16505.

84. Walsh SK, Imrie RM, Matuszewska M, Paterson GK, Weinert LA, Hadfield JD, et al. The host phylogeny determines viral infectivity and replication across Staphylococcus host species. PLoS Pathog. 2023;19(6):e1011433. https://doi.org/10.1371/journal.ppat.1011433 PMID: 37289828

85. Holtappels D, Alfenas-Zerbini P, Koskella B. Drivers and consequences of bacteriophage host range. FEMS Microbiol Rev. 2023;47(4):fuad038.
https://doi.org/10.1093/femsre/fuad038 PMID: 37422441

86. Alseth EO, Roush C, Irby I, Kopylov M, Bobe D, Diggs MW, et al. Mystique, a broad host range Acinetobacter phage, reveals the impact of culturing conditions on phage isolation and infectivity. PLoS Pathog. 2025;21(4):e1012986. https://doi.org/10.1371/journal.ppat.1012986 PMID: 40208916

87. Sprouffske K, Wagner A. Growthcurver: an R package for obtaining interpretable metrics from microbial growth curves. BMC Bioinformatics. 2016;17:172. https://doi.org/10.1186/s12859-016-1016-7 PMID: 27094401

88. Blazanin M, Vasen E, Vilaró Jolis C, An W, Turner PE. Quantifying phage infectivity from characteristics of bacterial population dynamics. Proc Natl Acad Sci U S A. 2025;122(37):e2513377122. https://doi.org/10.1073/pnas.2513377122 PMID: 40920928

89. Almeida-Neto M, Guimarães P, Guimarães PR Jr, Loyola RD, Ulrich W. A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement. Oikos. 2008;117(8):1227-39. https://doi.org/10.1111/j.0030-1299.2008.16644.x
90. Galeano J, Pastor JM, Iriondo JM. Weighted-interaction nestedness estimator (WINE): A new estimator to calculate over frequency matrices. Environ Model Softw. 2009;24(11):1342-6.

91. Newman MEJ, Girvan M. Finding and evaluating community structure in networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2004;69(2 Pt 2):026113. https://doi.org/10.1103/PhysRevE.69.026113 PMID: 14995526

92. Reichardt J, Bornholdt S. Statistical mechanics of community detection. Phys Rev E. 2006;74(1):016110.

93. Traag VA, Bruggeman J. Community detection in networks with positive and negative links. Phys Rev E. 2009;80(3):036115.

94. Clauset A, Newman MEJ, Moore C. Finding community structure in very large networks. Phys Rev E. 2004;70(6):066111.

95. Newman MEJ. Finding community structure in networks using the eigenvectors of matrices. Phys Rev E. 2006;74(3):036104.

96. Beckett SJ. Improved community detection in weighted bipartite networks. R Soc Open Sci. 2016;3(1):140536. https://doi.org/10.1098/rsos.140536 PMID: 26909160

97. Blondel VD, Guillaume J-L, Lambiotte R, Lefebvre E. Fast unfolding of communities in large networks. J Stat Mech. 2008;2008(10):P10008. https://doi.org/10.1088/1742-5468/2008/10/p10008

98. Rodríguez-Gironés MA, Santamaría L. A new algorithm to calculate the nestedness temperature of presence-absence matrices. J Biogeogr. 2006;33(5):924-35.

99. Moury B, Fabre F, Hébrard E, Froissart R. Determinants of host species range in plant viruses. J Gen Virol. 2017;98(4):862-73. https://doi.org/10.1099/jgv.0.000742 PMID: 28475036

100. Weitz JS, Poisot T, Meyer JR, Flores CO, Valverde S, Sullivan MB, et al. Phage-bacteria infection networks. Trends Microbiol. 2013;21(2):82-91. https://doi.org/10.1016/j.tim.2012.11.003 PMID: 23245704

101. Anderson MJ. Permutational multivariate analysis of variance (PERMANOVA). In: Wiley StatsRef: Statistics Reference Online [Internet]. John Wiley & Sons, Ltd; 2017 [cited 2025 June 4]. p. 1-15. Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/9781118445112.stat07841

102. Dixon P. Vegan, a package of R functions for community ecology. J Veg Sci. 2003;14(6):927-30.

103. Arbizu PM. pmartinezarbizu/pairwiseAdonis [Internet]. 2025 [cited 2025 Apr 14]. Available from: https://github.com/pmartinezarbizu/pairwiseAdonis

104. Brooks ME, Kristensen K, Benthem KJ, van Magnusson A, Berg CW, et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 2017;9(2):378.

105. Allen B, Kon M, Bar-Yam Y. A new phylogenetic diversity measure generalizing the shannon index and its application to phyllostomid bats. Am Nat. 2009;174(2):236-43. https://doi.org/10.1086/600101 PMID: 19548837