Genomic and in vitro characterization of two lytic bacteriophages infecting multidrug-resistant Erwinia sp. strain AnSW2-5

TL;DR

Two new bacteriophages were studied for their ability to kill drug-resistant Erwinia bacteria, showing strong effectiveness and resistance suppression.

Contribution

Discovery and characterization of two lytic bacteriophages with synergistic effects against MDR Erwinia sp.

Findings

Phage cocktail reduced bacteria by >80% and suppressed growth for 72 hours.

Dual-phage treatment prevented resistant mutant emergence with FoR < 10⁻⁸.

P-A and P-K phages showed distinct lytic properties, including latent period and burst size differences.

Abstract

Multidrug-resistant (MDR) Erwinia pathogens present a critical threat to global crop health and agricultural sustainability. In this study, we characterized two novel lytic bacteriophages, AnSW2-5-P-A (family Autographiviridae) and AnSW2-5-P-K (class Caudoviricetes), targeting the MDR Erwinia sp. strain AnSW2-5. Comparative genomics and TEM analysis revealed distinct virion architectures and confirmed the absence of lysogeny-associated genes, ensuring their safety as biocontrol agents. One-step growth assays demonstrated that P-A has a shorter latent period (~ 20 min), while P-K exhibits a significantly larger burst size (~ 110 PFU/cell). In co-culture assays, the dual-phage cocktail demonstrated a profound synergistic effect, achieving > 80% bacterial reduction (p < 0.05) and maintaining sustained suppression of the host for 72 h. Notably, the cocktail effectively prevented the emergence of resistant mutants, reducing the frequency of resistance (FoR) to below the limit of detection (< 10⁻⁸). These findings highlight the potential of using genetically diverse phage pairs with complementary lytic activities as a robust, resistance-suppressive biocontrol strategy against MDR phytopathogens.

The online version contains supplementary material available at 10.1038/s41598-026-39563-8.

Figures

Phylogenetic tree based on 16 S rRNA gene sequences showing the relationship of strain *Erwinia* sp. AnSW2-5 to closely related taxa within the genus *Erwinia* and related genera including *Pantoea*, *Winslowiella*, and *Phytobacter*. The tree was constructed using both neighbor-joining (NJ) and maximum likelihood (ML) methods based on the Jukes–Cantor model with 1,000 bootstrap replicates. Bootstrap values ≥ 70% (NJ/ML) are shown at the nodes. *Budvicia aquatica* DSM 5075ᵀ was used as the outgroup. Strain AnSW2-5 clustered most closely with *Erwinia rhapontici* ATCC 29283ᵀ with high bootstrap

Isolation and morphological characterization of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K. (A) Schematic diagram of the phage enrichment and isolation process from freshwater samples using a double-layer agar method. (B) Spot assay showing clear lysis zones on *Erwinia* sp. AnSW2-5 lawn. (C) Transmission electron micrograph (TEM) of AnSW2-5-P-A, exhibiting an icosahedral capsid and short tail typical of *Autographiviridae*. (D) TEM image of AnSW2-5-P-K, showing an icosahedral head and long non-contractile tail consistent with *Caudoviricetes* (*Myoviridae*-like) morphology.

Circular genome maps of Erwinia phages AnSW2-5-P-A (A) and AnSW2-5-P-K (B), generated using CGView. CDS are shown in blue. The inner rings depict GC content and GC skew. Annotated genes include structural proteins, DNA packaging elements, lytic enzymes, and host interaction factors. The modular genomic organization reflects typical *Caudoviricetes* phages and supports their lytic nature.

One-step growth curves of phages AnSW2-5-P-A and AnSW2-5-P-K. The infection kinetics of the two phages were determined using *Erwinia* sp.

Synergistic antibacterial activity and suppression of resistance emergence. (A, B) Bacterial growth curves of *Erwinia* sp. AnSW2-5 treated with single phages (P-A, P-K) and their cocktail over 72 h. (C) Frequency of resistance (*FoR*) to single phages and the phage cocktail. The red dashed line represents the limit of detection (LOD, 10^− 8^). The dual-phage cocktail maintained complete growth inhibition by significantly reducing the emergence of resistant variants. Error bars indicate SD from three independent replicates.

Full text

Introduction

The increasing prevalence of antibiotic-resistant bacteria presents a critical challenge to global agriculture, threatening crop productivity, food security, and public health alike^1–4^. In particular, the extensive and often unselective use of antibiotics in farming systems has accelerated the emergence of multidrug-resistant (MDR) pathogens, which has led to treatment failures and recurring outbreaks in plant systems^5–7^. Notably, members of the Erwinia genus, including Erwinia amylovora and Erwinia carotovora, are among the most destructive phytopathogens, causing fire blight and soft rot in economically valuable crops such as apples, pears, and potatoes^8–12^.

Current management strategies rely heavily on antibiotics (e.g., streptomycin) and copper-based compounds. However, these methods have become increasingly ineffective due to the rapid evolution of resistant strains^13–15^. For instance, streptomycin-resistant E. amylovora has been reported with rising frequency across orchards in North America and Europe, significantly reducing the efficacy of conventional disease control programs^11,12,16^. This highlights the need for sustainable biocontrol methods effective against MDR phytopathogens and compatible with eco-friendly practices^17^.

Bacteriophages (phages), viruses that specifically infect and lyse bacterial cells, are gaining attention as promising biocontrol agents for MDR bacteria^18–20^. Owing to their high specificity, phages can selectively target pathogens without disrupting beneficial microbiota^21,22^. Additionally, their evolutionary adaptability provides a dynamic advantage over conventional antimicrobials, which tend to exert static selective pressure and promote resistance^23,24^. While phage therapy has been explored in human and veterinary medicine, its application in plant disease management remains limited^25,26^. However, the commercial availability of EPA-registered phage biopesticides such as AgriPhage^®^ for fire blight and bacterial spot control in apples, pears, and tomatoes demonstrates growing practical interest in agricultural phage applications^27,28^.

To date, several phages infecting Erwinia species have been isolated. However, most exhibit narrow host ranges, reduced efficacy under environmental stresses such as UV exposure or desiccation, and limited persistence in field conditions^29–31^. Moreover, challenges related to phage stability, bacterial resistance evolution, and delivery strategy optimization continue to impede their broad implementation in plant disease control^31^. Recent studies suggest that phage cocktails comprising genetically distinct bacteriophages can broaden host range and enhance bactericidal efficacy. In particular, a combination of phages AnSW2-5-P-A and AnSW2-5-P-K showed greater suppression of MDR Erwinia than single-phage treatments, indicating complementary and synergistic interactions^32^. This approach is further supported by a recent European study, where a phage cocktail demonstrated effective control of E. amylovora in planta, reinforcing the feasibility of phage-based biocontrol under realistic agricultural conditions^33^.

In this study, we report the isolation and comprehensive in vitro characterization of two novel lytic bacteriophages, AnSW2-5-P-A and AnSW2-5-P-K, that specifically target the MDR Erwinia sp. strain AnSW2-5. We examined their morphological and genomic features, infection kinetics (including one-step growth analysis), and bacteriolytic activities. We hypothesized that the genetic divergence between these two phages—belonging to the Autographiviridae family and the class Caudoviricetes—enables complementary antibacterial effects. Our findings provide a genomic and experimental basis for controlling MDR plant pathogens, focusing on the potential of specific phage pairs to prevent resistance emergence in vitro.

Materials and methods

Isolation and 16 S rRNA gene sequencing of host bacteria

Soil samples were collected on 26 April 2018 from a cattle kraal (37°54′30.3″N, 128°48′27.1″E), a site previously identified in our field survey as being frequently impacted by livestock waste and antibiotic use, and therefore considered a hotspot for MDR bacteria^34^. The collected samples were surface-sterilized using 70% ethanol followed by rinsing with sterile distilled water, then homogenized and plated onto Reasoner’s 2 A (R2A) agar (MBcell, Seoul, South Korea). No plant materials were collected.

For selective enrichment, R2A agar was supplemented with penicillin–streptomycin (100 U/mL–100 µg/mL), tetracycline (10 µg/mL), lincomycin (10 µg/mL), and novobiocin (25 µg/mL) (Sigma-Aldrich, St. Louis, MO, USA). Antibiotic concentrations were chosen based on prior studies employing these drugs for the selective isolation of environmental bacteria under antibiotic pressure^35^. Plates were incubated at 20 °C for up to 14 days. Colony emergence was delayed under low-nutrient conditions with antibiotic selection, with initial colonies typically appearing after 5–7 days. Therefore, daily monitoring was extended to a maximum of 14 days to ensure recovery of slow-growing or stressed cells. Well-isolated colonies that grew on antibiotic-containing media were purified and designated Erwinia sp. strain AnSW2-5.

Antibiotic susceptibility of strain AnSW2-5 was tested by the disk diffusion method on Mueller–Hinton agar (Sigma-Aldrich). Commercial disks contained ampicillin 10 µg, cefotaxime 30 µg, penicillin 10 IU, vancomycin 30 µg, erythromycin 15 µg, kanamycin 30 µg, streptomycin 10 µg, tetracycline 30 µg, rifampin 5 µg, sulfamethoxazole 23.75 µg, trimethoprim 5 µg, chloramphenicol 30 µg, gentamicin 10 µg, colistin 10 µg, novobiocin 30 µg, and nalidixic acid 30 µg. Because no CLSI or EUCAST interpretive criteria exist for phytopathogenic Erwinia spp., inhibition-zone diameters are reported as raw values and interpreted descriptively rather than by applying clinical breakpoints. The strain was stored at − 80 °C in 20% (v/v) glycerol and deposited in the Freshwater Bioresources Culture Collection (FBCC) under accession number FBCC-B5912.

Genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen). The nearly complete 16 S rRNA gene was amplified using universal primers (27 F, 518 F, 800R, 1492R) under standard PCR cycling conditions^36^. Taxonomic assignment was performed using the EzBioCloud database^37^, which provides a curated and regularly updated repository of type-strain sequences with integrated OrthoANI metrics, thereby reducing misidentification risk. Phylogenetic trees were constructed using neighbor-joining (NJ)^38^ and maximum-likelihood (ML)^39^ methods in MEGA version 7.0^40^ with the Jukes–Cantor model; bootstrap analysis was performed with 1,000 replicates. The 16 S rRNA gene of strain AnSW2-5 showed 99.7% similarity to Erwinia rhapontici ATCC 29,283, but the OrthoANI value to that species was 91.6%, below the accepted threshold (≥ 95–96%)^41^. Thus, strain AnSW2-5 was designated as Erwinia sp., most closely related to E. rhapontici.

Genome sequencing, assembly, and annotation of host bacteria

Genomic DNA of strain AnSW2-5 was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. DNA quality and purity were assessed using a NanoDrop spectrophotometer by measuring the A260/A280 and A260/A230 ratios, and integrity was confirmed by 1% agarose gel electrophoresis.

Whole-genome sequencing was conducted using the Illumina HiSeq 2000 platform (DNA Link, Seoul, South Korea), generating high-throughput paired-end reads. Raw reads were quality-trimmed and filtered using Trimmomatic v0.39 to remove adapter sequences and low-quality bases^42^. De novo assembly was performed using CLC Genomics Workbench v9.0 (Qiagen) under stringent parameters: minimum read length of 100 bp, Phred quality score cutoff of Q30, and expected genome coverage of > 100×.

Genome annotation was carried out using the Rapid Annotation using Subsystem Technology (RAST) server^43^. Predicted coding sequences (CDSs) were functionally validated through BLASTP searches against the NCBI non-redundant (nr) protein database, using an E-value threshold of 1e^5^ and a minimum sequence identity of 30%, which is a standard threshold for detecting remote homology in sequence analysis. Functional categorization was further refined through pathway mapping using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database via the KEGG Mapper tool^44^.

To identify antibiotic resistance genes (ARGs), the assembled genome was screened using ResFinder 4.0^45^, the Comprehensive Antibiotic Resistance Database (CARD)^46^, and AMRFinderPlus^47^. ResFinder was used to detect acquired resistance genes, applying thresholds of ≥ 90% sequence identity and ≥ 60% coverage. CARD analysis was performed with the Resistance Gene Identifier (RGI) tool to detect both horizontally transferred resistance genes and chromosomal mutations. AMRFinderPlus provided complementary predictions based on both nucleotide and amino acid sequences. Cross-validation across these three tools was performed to ensure consistency and accuracy.

Phage isolation, purification, and concentration

Bacteriophages targeting Erwinia sp. strain AnSW2-5 were isolated from freshwater samples collected at the same site as the host. For phage enrichment, 400 mL of freshwater filtrate (0.22 μm pore-size, Millipore, Burlington, MA, USA) was mixed with 100 mL of 5× R2A broth (BD Difco, Franklin Lakes, NJ, USA) and 30 mL of an exponentially growing host culture (OD_600_ ≈ 0.5)^48^. The mixture was incubated at 30 °C for seven days. To monitor phage propagation, 10 mL aliquots were collected at three-day intervals, treated with 2 mL of chloroform to lyse remaining bacteria, and centrifuged at 10,000 × g for 10 min to obtain clarified lysates^49^.

Phage purification was conducted via the double-layer agar overlay method. Serial 10-fold dilutions of the lysates were mixed with 0.5 mL of host culture and 6 mL of molten R2A top agar (0.7% w/v agar), then overlaid onto R2A bottom agar (1.5% w/v agar). Plates were incubated at 30 °C for up to one week; this extended incubation was essential to maximize the recovery of diverse phages, including those with slow plaque development due to longer latent periods or the slow growth of the host on R2A medium. Individual plaques were picked and resuspended in SM buffer (50 mM Tris-HCl pH 7.5, 100 mM NaCl, 8 mM MgSO_4_). This purification process was repeated for five successive rounds to ensure monoclonality. Based on distinct plaque morphologies, two phages were selected and designated AnSW2-5-P-A and AnSW2-5-P-K.

For high-purity preparations, phages were amplified in large-scale cultures (500 mL) until complete lysis occurred. Lysates were filtered (0.22 μm) and concentrated using 10% (w/v) polyethylene glycol (PEG)−8000 and 1 M NaCl, followed by overnight incubation at 4 °C. After centrifugation at 12,000 × g for 30 min, pellets were resuspended in SM buffer. Final purification was achieved by cesium chloride (CsCl) step-gradient ultracentrifugation using density layers of 1.3, 1.5, and 1.7 g/mL. The gradient was centrifuged at 100,000 × g for 4 h at 4 °C. The resulting opalescent phage bands were collected, dialyzed against SM buffer to remove CsCl, and stored for further analysis.

Phage morphological analysis

The structural morphology of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K was examined using transmission electron microscopy (TEM) to determine capsid architecture, tail structure, and viral family classification. Phage lysates were first purified by polyethylene glycol (PEG)-NaCl precipitation followed by ultracentrifugation at 30,000 × g for 1 h at 4 °C. For TEM analysis, purified phage suspensions were adsorbed onto 200-mesh carbon-coated copper grids (Electron Microscopy Sciences), negatively stained with 2% uranyl acetate for 10 s, and air-dried^50^. Phage morphology was visualized using a Philips CM200 TEM operating at 80 kV. Morphometric parameters, including capsid diameter and tail length, were measured for at least 10 randomly selected virions using ImageJ software. Mean values and standard deviations were calculated.

This morphological analysis provided key insights into virion symmetry, size, and tail architecture, enabling preliminary taxonomic classification according to the latest guidelines from the International Committee on Taxonomy of Viruses (ICTV)^51^.

Phage genome extraction, sequencing, and analysis

Phage genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following a modified protocol from Molecular Cloning: A Laboratory Manual^48^. Briefly, 70 µL of concentrated phage lysate was mixed with 30 µL of proteinase K, 3 µL of RNase A, and 300 µL of ATL buffer, followed by incubation at 56 °C to digest capsid proteins. Subsequently, 300 µL of AL buffer and 100% ethanol were added, and the mixture was loaded onto a Qiagen spin column. The column was washed with AW1 and AW2 buffers, and DNA was eluted using EB buffer. DNA quality and concentration were measured using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and integrity was verified by 1% agarose gel electrophoresis.

Whole-genome sequencing of phages AnSW2-5-P-A and AnSW2-5-P-K was conducted using the Illumina HiSeq 2500 platform (DNA Link, Seoul, South Korea), producing 150-bp paired-end reads. Raw reads were trimmed and filtered using Trimmomatic v0.39 to remove adapters and low-quality sequences^42^. High-quality reads were assembled de novo using SPAdes v3.15.3^52^, and assembly quality was assessed with QUAST^53^.

To ensure complete genome coverage, PCR-based gap filling was performed using primers designed with Primer3 targeting regions with < 10× coverage. PCR products were validated by Sanger sequencing under standard thermal cycling conditions (30 cycles, 55 °C annealing).

Open reading frames (ORFs) were predicted using a combination of RAST^43^, GeneMark v3.25^54^, and GLIMMER v3.02^55^, followed by manual curation with RAST and Artemis^56^. Functional annotation was performed via BLASTp against the NCBI non-redundant (nr) database, and supplemented by domain prediction using Pfam^57^, TIGRFAM^58^, and COG^59^. Additional tools included InterProScan^60^ for conserved domain identification, TMHMM^61^ for transmembrane helices, and SignalP v6.0^62^ for signal peptide prediction.

Phylogenetic analysis was conducted based on two conserved viral marker genes: the terminase large subunit (TerL) and DNA polymerase A. Multiple sequence alignments were performed using CLUSTALW, and phylogenetic trees were reconstructed using the maximum likelihood (ML) method in MEGA v7.0^40^, with the JTT matrix model and gamma distribution (Γ) for rate variation. Node robustness was assessed by 1,000 bootstrap replicates, and values ≥ 70% were considered statistically supported.

Phage one-step growth curve

To determine the latent period and burst size of phages AnSW2-5-P-A and AnSW2-5-P-K, a one-step growth experiment was performed as previously described^63^ with minor modifications. Exponentially growing Erwinia sp. AnSW2-5 (1 × 10^8^ CFU/mL) was infected with each phage at a multiplicity of infection (MOI) of 0.01. The mixture was incubated at 25 °C for 10 min to allow for phage adsorption, followed by centrifugation at 13,000 × g for 3 min to remove unadsorbed phages. The pellet was resuspended in 50 mL of fresh R2A broth and incubated at 30 °C. Samples were collected at 5–10 min intervals for up to 100 min. Phage titers were determined using the double-layer agar method. The latent period was defined as the time interval between adsorption and the beginning of the first rise in phage titer. The burst size was calculated as the ratio of the final count of liberated phage particles to the initial count of infected bacterial cells.

Analysis of phage proliferation and bacteriolytic activity in co-culture

The proliferation dynamics and lytic activity of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K against Erwinia sp. AnSW2-5 were evaluated through co-culture experiments. The multiplicity of infection (MOI) was set at 1.0 to ensure an equal initial ratio of phages to bacterial cells. The phage–host mixture was incubated at 25 °C for 10 min to facilitate adsorption, followed by centrifugation at 13,000 × g for 3 min to remove unadsorbed phages. The resulting pellet was resuspended in R2A broth and incubated at 30 °C to initiate co-culture (defined as time point 0 h).

Bacterial lysis was monitored by measuring optical density at 600 nm (OD_600_) every 2 h for a total of 72 h using a Novaspec Pro spectrophotometer (Biochrom, UK). In parallel, bacterial viability and phage proliferation were evaluated at each time point by quantifying colony-forming units (CFU/mL) and plaque-forming units (PFU/mL), respectively.

To compare the efficacy of different treatments, three infection conditions were tested: (i) AnSW2-5-P-A alone, (ii) AnSW2-5-P-K alone, and (iii) a combination of both phages. A two-way analysis of variance (ANOVA) was conducted with phage treatment and time as independent variables, and bacterial reduction (measured via OD_600_ and CFU/mL) as the dependent variable. Tukey’s post hoc test (p < 0.05) was applied for multiple comparisons among groups. All experiments were conducted in triplicate with three independent biological replicates.

Determination of frequency of resistance (FoR)

The frequency of bacteriophage-resistant mutants was determined according to previously described methods^64^ with slight modifications. Briefly, a bacterial culture of Erwinia sp. AnSW2-5 was grown to the early exponential phase (≈ 10^8^ CFU/mL). For single-phage treatments (AnSW2-5-P-A or AnSW2-5-P-K) and the dual-phage cocktail, phages were mixed with the bacterial suspension at a multiplicity of infection (MOI) of 10. The mixtures were incubated at 25 °C for 15 min to allow phage adsorption and then serially diluted in sterile PBS. The dilutions were plated onto R2A agar plates using the double-layer agar method and incubated at 25 °C for 48 h. The FoR was calculated by dividing the number of resistant colonies appearing on the phage-treated plates by the total number of viable bacteria in the non-phage-treated control group. The limit of detection (LOD) was 10^− 8^. All experiments were performed in triplicate.

Results

Isolation and 16 S rRNA gene phylogeny of host bacteria

Soil samples collected on April 26, 2018, from a cattle kraal near the Nakdong River, South Korea, were processed for bacterial isolation. After surface sterilization and plating on R2A agar with and without antibiotics, colonies capable of growing under selective pressure were isolated and screened for multidrug resistance.

The isolate AnSW2-5 exhibited resistance to multiple antibiotics. Disk diffusion assays showed inhibition zones of 5.2 ± 0.5 mm (cefotaxime), 7.1 ± 0.4 mm (ampicillin), 6.8 ± 0.3 mm (penicillin), and complete resistance to vancomycin (no zone of inhibition). Furthermore, the strain showed no inhibition zones (0.0 ± 0.0 mm) for streptomycin, kanamycin, and tetracycline, indicating high-level resistance to these classes. Detailed susceptibility profiles for 13 tested antibiotics and their corresponding resistance determinants are summarized in Table 1. These findings confirm a strong MDR phenotype, consistent with previous reports on antibiotic-resistant Erwinia species.

Table 1. Antibiotic susceptibility profile and corresponding resistance genes of Erwinia sp. AnSW2-5.Antibiotic ClassAntibiotic (Disk content)Inhibition Zone Diameter (mm)Identified Resistance Determinantsβ-lactamsAmpicillin (10 µg)7.1 ± 0.4bla_TEM−1_, acrAB-tolCPenicillin (10 IU)6.8 ± 0.3bla_TEM−1_Cefotaxime (30 µg)5.2 ± 0.5bla_TEM−1_AminoglycosidesStreptomycin (10 µg)0.0 ± 0.0 (R)aph(3’)-Ia,* aac(6’)-Ib3Kanamycin (30 µg)0.0 ± 0.0 (R)aph(3’)-Ia, aac(6’)-Ib3Gentamicin (10 µg)8.5 ± 0.6aac(6’)-Ib3*, acrAB-tolCTetracyclinesTetracycline (30 µg)0.0 ± 0.0 (R)tetA,* tetR*,* acrAB-tolC*PhenicolsChloramphenicol (30 µg)9.2 ± 0.4 acrAB-tolC QuinolonesNalidixic acid (30 µg)10.5 ± 0.7 acrAB-tolC GlycopeptidesVancomycin (30 µg)0.0 ± 0.0 (R)(Intrinsic resistance)OthersRifampin (5 µg)7.8 ± 0.5 acrAB-tolC Sulfamethoxazole (23.75 µg)12.4 ± 0.8-Novobiocin (30 µg)11.2 ± 0.6-

The nearly full-length 16 S rRNA gene sequence of AnSW2-5 was obtained and compared to reference sequences in the GenBank and EzBioCloud databases using BLAST. The sequence exhibited 99.7% identity to Erwinia rhapontici ATCC 29,283. Phylogenetic trees were constructed using both neighbor-joining (NJ) and maximum likelihood (ML) methods in MEGA v7.0, with evolutionary distances calculated by the Jukes–Cantor model. Both analyses strongly supported the close relationship between AnSW2-5 and E. rhapontici, with bootstrap values of 99% (NJ) and 98% (ML), confirming its taxonomic assignment within the genus Erwinia (Fig. 1).

Fig. 1. Phylogenetic tree based on 16 S rRNA gene sequences showing the relationship of strain Erwinia sp. AnSW2-5 to closely related taxa within the genus Erwinia and related genera including Pantoea, Winslowiella, and Phytobacter. The tree was constructed using both neighbor-joining (NJ) and maximum likelihood (ML) methods based on the Jukes–Cantor model with 1,000 bootstrap replicates. Bootstrap values ≥ 70% (NJ/ML) are shown at the nodes. Budvicia aquatica DSM 5075ᵀ was used as the outgroup. Strain AnSW2-5 clustered most closely with Erwinia rhapontici ATCC 29283ᵀ with high bootstrap support (99% NJ/98% ML).

Genome analysis of Erwinia sp. AnSW2-5

The genome of Erwinia sp. AnSW2-5 was sequenced using the Illumina HiSeq 2000 platform, yielding a draft assembly consisting of 62 contigs with an average coverage depth of 115×. The total genome size was 5,032,288 bp with a G + C content of 54.5 mol%, which is consistent with reported values for other Erwinia genomes.

A total of 4,664 protein-coding sequences and 66 RNA genes were predicted. Functional categorization using the RAST annotation server assigned coding sequences to various subsystem categories, including carbohydrate metabolism, protein metabolism, stress response, and virulence (Figure S1). Subsystem coverage analysis indicated that 34% of annotated features were associated with known biological functions, while 66% remained unclassified.

To elucidate the genetic basis of the observed multidrug resistance phenotype, the genome was screened using ResFinder 4.1 and the Comprehensive Antibiotic Resistance Database (CARD). Genomic analysis of Erwinia sp. AnSW2-5 revealed multiple open reading frames (ORFs) encoding resistance determinants that correlate with the MDR phenotype. As summarized in Table 1, the identified genetic repertoire includes the β-lactamase gene bla_TEM-1_, aminoglycoside-modifying enzymes (aph(3′)-Ia and aac(6′)-Ib3), and tetracycline resistance determinants (tetA and tetR). These genomic features are consistent with the phenotypic resistance profiles observed in disk diffusion assays, where the strain exhibited significant resistance to the corresponding antibiotic classes. In particular, the identification of the acrAB-tolC system—a tripartite multidrug efflux pump—suggests a robust mechanism for broad-spectrum resistance by actively extruding diverse antimicrobial agents^65^. The co-occurrence of these site-specific modifying enzymes and generalized efflux systems underscores the complex resistance strategy of strain AnSW2-5, providing a genetic explanation for its survival under the high antibiotic pressure of its source environment.

Genomic relatedness was assessed using Orthologous Average Nucleotide Identity (OrthoANI) analysis against five type strains of Erwinia species. The highest ANI value was observed with E. rhapontici ATCC 29283ᵀ (91.6%), followed by E. persicina NBRC 102418ᵀ (83.8%), E. aphidicola USMM130 (81.8%), E. tasmaniensis Et1/99ᵀ (80.7%), and E. billingiae Eb661ᵀ (77.8%) (Figure S2). All ANI values were below the 95–96% species-level threshold. According to the taxonomic criteria proposed for prokaryotic species demarcation^66^, these results indicate that strain AnSW2-5 is genetically distinct from previously described Erwinia species. However, following the standard polyphasic approach, the strain was designated as Erwinia sp. AnSW2-5 until further phenotypic and chemotaxonomic characterization is completed.

Further genome mining revealed genes associated with environmental adaptation, including nitrogen metabolism (nif operon, urease cluster), oxidative stress resistance (katE, sodA), and biofilm formation (bcsA, pga operon)^67^. These genes likely support survival in soil and plant-associated niches by facilitating nitrogen assimilation, stress mitigation, and surface adherence.

Collectively, these findings demonstrate that Erwinia sp. AnSW2-5: (i) possesses extensive antibiotic resistance, (ii) is genomically distinct from known Erwinia species, and (iii) harbors genetic features associated with environmental resilience. These characteristics support its relevance as a model organism for bacteriophage-host interaction studies and its potential as a target for phage-based biocontrol strategies in agriculture.

Bacteriophage isolation and morphological characteristics

Bacteriophages infecting Erwinia sp. AnSW2-5 were isolated from freshwater samples collected at the same site as the host strain (Fig. 2A). Following enrichment and plaque purification, two distinct phages were identified and designated AnSW2-5-P-A and AnSW2-5-P-K (Fig. 2B).

Fig. 2. Isolation and morphological characterization of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K. (A) Schematic diagram of the phage enrichment and isolation process from freshwater samples using a double-layer agar method. (B) Spot assay showing clear lysis zones on Erwinia sp. AnSW2-5 lawn. (C) Transmission electron micrograph (TEM) of AnSW2-5-P-A, exhibiting an icosahedral capsid and short tail typical of Autographiviridae. (D) TEM image of AnSW2-5-P-K, showing an icosahedral head and long non-contractile tail consistent with Caudoviricetes (Myoviridae-like) morphology.

Transmission electron microscopy (TEM) revealed that AnSW2-5-P-A possesses an icosahedral capsid (55–65 nm) with a short, non-contractile tail (5–15 nm) (Fig. 2C, Figure S3). This morphology is characteristic of the family Autographiviridae^51,68^, previously grouped within the Podoviridae family^51,69^, known for rapid, single-step genome injection^70^.

In contrast, AnSW2-5-P-K exhibited a straight, non-flexible, and contractile tail (75–95 nm) with a distinct tail sheath and baseplate structure (Fig. 2D, Figure S4). The icosahedral capsid measured 60–70 nm. This structural organization aligns with the Myoviridae-like morphology within the class Caudoviricetes^51^. Members of this lineage are known to infect diverse gammaproteobacterial hosts, including Erwinia, utilizing a sophisticated tail contraction mechanism^71^. Detailed genome re-annotation (Table S2) supported this classification by identifying essential structural genes, including those encoding tail sheath, tube, and baseplate assembly proteins^72^. These findings collectively confirm that AnSW2-5-P-K possesses a complex, contractile infection apparatus typical of the Caudoviricetes.

The morphological differences between the two phages suggest distinct infection strategies. AnSW2-5-P-A, with its compact structure, likely facilitates rapid adsorption and genome injection, consistent with the behavior of Autographiviridae. In contrast, AnSW2-5-P-K’s contractile tail apparatus supports a more forceful and regulated genome delivery through the host cell envelope, reflecting traits typical of contractile-tailed phages in the Caudoviricetes.

The architectural divergence between the two phages may reflect distinct ecological adaptations. Compact, short-tailed phages such as AnSW2-5-P-A may be better suited to high-density, surface-associated microbial environments where rapid replication is advantageous. In contrast, complex, contractile-tailed phages like AnSW2-5-P-K may confer greater ecological versatility, potentially enhancing infection efficiency in diverse or spatially structured habitats, such as biofilms, where overcoming complex cell wall barriers is critical. Collectively, these findings indicate that AnSW2-5-P-A and AnSW2-5-P-K represent morphologically and functionally distinct bacteriophage lineages, likely shaped by divergent evolutionary and environmental selection pressures.

Genome analysis of Erwinia phages

The complete genomes of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K were sequenced using the Illumina HiSeq 2500 platform. AnSW2-5-P-A has a linear double-stranded DNA (dsDNA) genome of 17,911 bp with a G + C content of 44.0%, while AnSW2-5-P-K comprises a 93,348 bp linear dsDNA genome with a G + C content of 50.0%. The assembly quality was validated using QUAST, which indicated high-quality assemblies with average sequencing depths of 115× and 98× for AnSW2-5-P-A and AnSW2-5-P-K, respectively.

A total of 35 ORFs were predicted in the AnSW2-5-P-A genome (Table S1), among which 21 were functionally annotated via BLASTp and conserved domain searches. Annotated functions included genes encoding a terminase large subunit (TerL, ORF4), endolysin (ORF6), tail fiber (ORF9), tail collar (ORF10), and structural components such as a phospholipase and baseplate protein (Fig. 3A). These features are consistent with the genomic organization of Autographiviridae phages, supporting rapid lytic infection strategies^51,68^.

Fig. 3. Circular genome maps of Erwinia phages AnSW2-5-P-A (A) and AnSW2-5-P-K (B), generated using CGView. CDS are shown in blue. The inner rings depict GC content and GC skew. Annotated genes include structural proteins, DNA packaging elements, lytic enzymes, and host interaction factors. The modular genomic organization reflects typical Caudoviricetes phages and supports their lytic nature.

The genome of AnSW2-5-P-K encodes 152 ORFs (Table S2), showing a more complex modular architecture (Fig. 3B). Functional annotation revealed genes related to DNA replication (DNA polymerase A, ORF117), transcription (DNA-dependent RNA polymerase, ORF128), DNA repair (RecA, ORF116; DNA ligase OB, ORF111), and host regulation (TetR family transcriptional regulator, ORF42). Lytic and structural genes include tail sheath proteins (ORF 136), tail tube proteins (ORF 135), baseplate assembly components (ORF 7), and multiple tail fiber proteins (e.g., ORFs 11, 140). These genomic features provide definitive evidence for the contractile infection apparatus observed in TEM analysis^71^. Conserved domain analysis using Pfam, TIGRFAM, and COG databases confirmed the presence of key functional modules, including DNA packaging (TerL, ORF58), structural proteins (ORFs 133, 135), and tail components (SLT domain, Gp138 N-terminal). Signal peptides and membrane-spanning regions were also predicted, such as in ORFs 55, 82, and 92, which may contribute to host interaction or viral release.

Phylogenetic relationships were inferred using ML trees based on the TerL gene (AnSW2-5-P-A) and DNA polymerase gene (AnSW2-5-P-K). AnSW2-5-P-A grouped with Kosakonia phage vB_KsaM-C1, forming a distinct lineage among Klebsiella and Pantoea phages (Figure S5). AnSW2-5-P-K clustered with Erwinia phages Zoomie, Loshitsaa2, and Micant with strong bootstrap support (99%), confirming its affiliation within Erwinia-specific Caudoviricetes phages (Figure S6).

Importantly, neither genome encodes integrase, recombinase, or repressor genes, indicating a strictly lytic lifestyle^73,74^. AnSW2-5-P-A, with a compact genome and T7-like architecture, is likely optimized for fast adsorption and genome delivery. In contrast, AnSW2-5-P-K exhibits a more complex genomic configuration utilizing a sophisticated tail contraction mechanism supported by the presence of tail sheath and baseplate genes (Table S2), potentially enabling broader ecological adaptability. These complementary genomic characteristics of AnSW2-5-P-A and AnSW2-5-P-K highlight their distinct yet synergistic biocontrol potential against multidrug-resistant Erwinia species.

One-step growth characteristics

To evaluate the lytic efficiency and developmental characteristics of the isolated phages, one-step growth experiments were performed using Erwinia sp. AnSW2-5 as the host (Fig. 4). The two phages exhibited distinct infection kinetics, reflecting their different taxonomic classifications and genomic architectures.

Fig. 4. One-step growth curves of phages AnSW2-5-P-A and AnSW2-5-P-K. The infection kinetics of the two phages were determined using Erwinia sp.

Phage AnSW2-5-P-A (Autographiviridae) displayed a rapid infection cycle characterized by a short latent period of approximately 20 min. The rise period continued until 40 min, reaching a plateau with a calculated burst size of approximately 70 ± 5 PFU per infected cell. This rapid generation time suggests a highly efficient host-takeover mechanism, likely facilitated by its streamlined genome and potent lytic enzymes^75^. In contrast, phage AnSW2-5-P-K (Caudoviricetes) exhibited a relatively prolonged latent period of approximately 35 min. However, it demonstrated a more robust proliferative capacity, with a rise period extending to 60 min and a significantly higher burst size of approximately 110 ± 8 PFU per infected cell compared to P-A^76^. This higher progeny yield aligns with the more complex genomic configuration and sophisticated infection apparatus identified in its contractile-tailed morphology.

These results indicate that while AnSW2-5-P-A provides a faster initial response to host availability, AnSW2-5-P-K contributes more substantially to the amplification of the phage population. This temporal staggering and the difference in progeny yield provide a mechanistic basis for the complementary lytic activity and sustained bacterial suppression observed in the subsequent co-culture assays^77^.

Bacteriolytic activity and synergistic effects in co-culture

The lytic activity of bacteriophages AnSW2-5-P-A and AnSW2-5-P-K was evaluated through co-culture experiments at a multiplicity of infection (MOI) of 1.0. In the untreated control group, Erwinia sp. AnSW2-5 exhibited steady growth, reaching a peak OD₆₀₀ of 2.7 at 72 h (Fig. 5A). In contrast, cultures treated with individual phages showed significant growth inhibition. AnSW2-5-P-A initiated detectable bacterial lysis around 32 h, whereas AnSW2-5-P-K exhibited a slightly delayed onset of inhibition. By 72 h, the OD_600_ in cultures treated with AnSW2-5-P-A and AnSW2-5-P-K was reduced by approximately 55% and 45%, respectively, relative to the control (Fig. 5B).

Fig. 5. Synergistic antibacterial activity and suppression of resistance emergence. (A,** B**) Bacterial growth curves of Erwinia sp. AnSW2-5 treated with single phages (P-A, P-K) and their cocktail over 72 h. (C) Frequency of resistance (FoR) to single phages and the phage cocktail. The red dashed line represents the limit of detection (LOD, 10^− 8^). The dual-phage cocktail maintained complete growth inhibition by significantly reducing the emergence of resistant variants. Error bars indicate SD from three independent replicates.

Notably, the dual-phage cocktail (AnSW2-5-P-A + AnSW2-5-P-K) demonstrated the most potent bacteriolytic effect, reducing the OD_600_ by over 80%. This synergistic interaction likely results from the targeting of non-overlapping host receptors, which broadens the overall host-killing spectrum and prevents the rapid emergence of resistance^78^. Statistical analysis using two-way ANOVA confirmed significant differences between the treatment groups (p < 0.05). Tukey’s post-hoc test further revealed that the dual-phage combination led to significantly greater bacterial suppression than either AnSW2-5-P-A (p < 0.01) or AnSW2-5-P-K (p < 0.05) alone. These findings were corroborated by viable cell count (CFU/mL) assays, which closely mirrored the OD_600_ results, a trend consistent with previous studies on phage-phage synergy against Erwinia species^79^. Concurrently, phage titers (PFU/mL) increased from 10^5^ to 10^9^ PFU/mL over the 72 h period, reflecting robust in situ phage replication.

To further elucidate the basis of the sustained growth inhibition, the frequency of resistance (FoR) was determined (Fig. 5C). In single-phage treatments, resistant variants emerged with FoR values of approximately 10^− 6^ and 10^− 7^ for P-A and P-K, respectively. These values are within the typical range reported for the emergence of spontaneous phage resistance in Gram-negative bacteria^80^. In contrast, the dual-phage cocktail markedly suppressed the emergence of resistant mutants, with FoR values falling below the limit of detection (< 10^− 8^). Consistent with these results, bacterial regrowth was observed in single-phage treatments, whereas the cocktail-treated group maintained complete inhibition throughout the 72 h period^78^. This synergistic effect indicates that the combination imposes a substantially higher genetic barrier to bacterial escape, as the probability of a single cell simultaneously acquiring multiple resistance mutations is drastically reduced^78,81^. These findings highlight the potential of phage cocktails for controlling multidrug-resistant Erwinia species.

Discussion

In this study, we characterized two novel lytic phages, AnSW2-5-P-A and AnSW2-5-P-K, exhibiting potent activity against an MDR Erwinia strain. Their divergent genomic architectures and infection strategies highlight their potential as complementary biocontrol agents in systems challenged by antibiotic resistance.

Phage AnSW2-5-P-A (Autographiviridae) employs a streamlined infection strategy with a short latent period (~ 20 min), likely facilitated by its functional endolysin^82^. In contrast, AnSW2-5-P-K (Caudoviricetes) possesses an expanded structural repertoire, including tail sheath and baseplate assembly components (ORFs 136 and 7). This contractile tail apparatus accounts for its prolonged latent period (~ 35 min) and significantly higher burst size (~ 110 PFU/cell), enabling sustained inhibitory effects^76^.

The dual-phage cocktail demonstrated a profound synergistic effect, maintaining complete growth inhibition for 72 h while single-phage treatments allowed rapid regrowth due to high resistance frequencies (10^− 6^ to 10^− 7^). Crucially, the cocktail suppressed the emergence of resistant variants to below the limit of detection (< 10^− 8^). This suggests that the combination of a short-tailed and a contractile-tailed phage imposes a higher genetic barrier to bacterial escape by creating non-overlapping infection bottlenecks^78,80^.

From a biosafety perspective, the absence of lysogeny-associated genes in both genomes is critical for agricultural application. While AnSW2-5-P-K’s genomic features (e.g., recA, DNA ligase) suggest robustness under environmental stress, we emphasize that further in vivo validation in plant-host models is essential to conclusively establish field efficacy. Future research should focus on protective formulations, such as microencapsulation, to translate this laboratory-scale success into reliable field-level biocontrol for phytopathogenic Erwinia species^83^.

Conclusions

This study demonstrates the potent biocontrol potential of two novel lytic bacteriophages, AnSW2-5-P-A and AnSW2-5-P-K, against a MDR Erwinia strain. Both phages exhibited complementary infection behaviors, characterized by rapid lysis by the Autographiviridae member (AnSW2-5-P-A) and a more regulated infection process by the Caudoviricetes member (AnSW2-5-P-K). Notably, the dual-phage cocktail not only enhanced bacterial suppression but also effectively suppressed the emergence of phage-resistant mutants to below the limit of detection (< 10^− 8^). These findings underscore the feasibility of phage-based biocontrol as a sustainable alternative to conventional chemical antibiotics in agricultural disease management, particularly for targeting MDR phytopathogens.

Despite these promising in vitro results, further research is essential to validate these findings in plant-host models (in vivo) to ensure long-term field efficacy and ecological safety. Key areas for future investigation include the impact of phage application on non-target microbial communities, the dynamics of phage resistance evolution, and the stability of phage formulations under fluctuating environmental conditions. Proactive monitoring of resistant bacterial variants and the strategic development of adaptive phage cocktails or sequential application schemes will be essential to maintain long-term control efficacy.

Future efforts should focus on optimizing formulation and delivery systems (e.g., microencapsulation), conducting comprehensive in-situ ecological impact assessments, and performing large-scale field trials across diverse agricultural landscapes. Addressing these practical challenges will be crucial for integrating phage-based biocontrol into mainstream integrated pest management (IPM) frameworks, ultimately advancing the goals of sustainable and environmentally responsible agriculture.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1