The genome and proteomes sequences of the analyzed phages were obtained from the NCBI Reference Sequence database: Tequintavirus T5 (GenBank, accession no. NC_005859.1), Lambdavirus lambda (GenBank, accession no. NC_001416.1), Dhillonvirus JLBYU60 (GenBank, accession no. NC_073052.1), Traversvirus tv24B (GenBank, accession no. NC_027984.1), Littlefixvirus 98PfluR60PP (GenBank, accession no. NC_070866.1).

Identity matrices for the analyzed proteins were generated using Clustal v2.1 (Madeira et al., 2019). Protein homology and sequence identity were assessed with BLASTp v2.9.0 (Camacho et al., 2009). Functional annotation and identification of conserved domains were performed using the NCBI Conserved Domain Database (accessed 2023–2024, https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) (Wang et al., 2023) and InterProScan v5.62–94.0 (Jones et al., 2014). Distant homologs of the analyzed proteins were identified using the Foldseek server (accessed 2023–2024, https://search.foldseek.com/) against comprehensive structural databases (van Kempen et al., 2024).

Protein homology modeling was conducted via the SWISS-MODEL server¹ (Waterhouse et al., 2018). The new receptor-binding proteins (RBPs) and protein oligomers were modeled using AlphaFold2 v2.3.2 (Jumper et al., 2021) and AlphaFold2-Multimer v3 (Evans et al., 2021), and implemented through ColabFold v1.5.2 (Evans et al., 2021; Jumper et al., 2021; Mirdita et al., 2022). The modeling protocol included template-based modeling with 30 recycles, and structure relaxation was performed using the AMBER (Assisted Model Building with Energy Refinement) option. The affinity of predicted protein structures to the planar outer membrane of Gram-negative bacteria was evaluated using the PPM 3.0 Web Server² (Lomize et al., 2022).

The Lambdavirus lambda protein NP_040600.1 was reconstructed by combining the N-terminal model—generated by SWISS-MODEL in its native homo-trimeric state—with a C-terminal model from AlphaFold2-Multimer v3. Assembly and refinement were carried out in ChimeraX v1.7 (Pettersen et al., 2021). Tail spike tip regions were digitally assembled using AlphaFold2-Multimer v3, restricting the input sequence to the final 200 amino acids of the C-terminus. Oligomeric states were inferred based on structural similarity to viral proteins with experimentally validated oligomeric conformations. The same procedure was applied to protein YP_010740310.1 from Dhillonvirus JLBYU60.

Protein-receptor interactions were predicted using heteromeric modeling in AlphaFold2-Multimer v3 with templates and 30 recycles. All multimeric predictions were manually curated to exclude biologically implausible interactions, particularly with membrane-embedded receptor regions inaccessible to RBPs. Modeling focused on phages with experimentally validated receptors: Dhillonvirus JLBYU60 and Traversvirus tv24B. For FhuA (UniProtKB ID: P06971), the full protein sequence was used, whereas for BamA (UniProtKB ID: P0A940), only the C-terminal β-barrel transmembrane domain (residues 427–810) was modelled.

The functional annotations were performed by Sphae snakemake workflow performed by Bhavya Papudeshi and colleagues (Papudeshi et al., 2025). This workflow allows made functional annotations by Pharokka (Bouras et al., 2023), Phold, which functional annotation using protein structures for determining the structural homology (unpublished), repository available³ and Phynteny, which is the synteny-based annotation of bacteriophage genes (Grigson et al., 2023). The protein sequences were also analysed by PHROG Sequence search for detection matching the Prokaryotic Virus Remote Homologous Groups (Terzian et al., 2021).

Predicted binding affinities of protein complexes were evaluated using the PRODIGY Web Server⁴ (Xue et al., 2016). Molecular dynamics (MD) simulations were conducted using GROMACS v2023.3 (Páll et al., 2020) with the CHARMM27 all-atom force field (Brooks et al., 1983, 2009). Simulations ran for 1 μs in a 0.15 M NaCl aqueous environment at 300 K. MD trajectories were analyzed using the MDAnalysis Python library v2.7.0 (Michaud-Agrawal et al., 2011), and root-mean-square deviation (RMSD) plots were generated with Plotly (Plotly Technologies Inc., Montréal, QC, 2015, https://plot.ly). Protein structures and cryo-EM maps were visualized and analyzed in ChimeraX v1.7, while MD trajectories were further examined using VMD v1.9 (Humphrey et al., 1996). Visualization videos were rendered using Adobe Premiere Elements 15 (Adobe Inc., San Jose, USA).

The bacteriophage Littlefixvirus 98PfluR60PP was propagated in a host bacterial culture. The culture medium was inoculated with the host strain and incubated at 25°C for approximately 3 h. Subsequently, the phage was added to the culture, followed by an additional 4-h incubation at 25°C to allow for host cell lysis. This procedure yielded a high-titer phage preparation, reaching approximately 1 × 10⁹ PFU/mL.

Following amplification, unlysed cells and cellular debris were removed by filtration. The resulting phage-containing supernatant was further purified and concentrated by cesium chloride (CsCl) ultracentrifugation at a final CsCl concentration of 1.5 g/mL. To prepare the gradient, 13.2 mL of phage suspension was gradually mixed with 9.9 g of CsCl (Thermo Scientific, J65950.22) and divided into three 4.4 mL ultracentrifuge tubes. Ultracentrifugation was performed for 23 h at 38,000 rpm (150,000 × g) at 4°C using a Thermo Scientific Sorvall™ WX + ultracentrifuge with the swinging bucket rotor (Thermo Scientific TH-660 rotor). The visible white-to-grey band containing phages was collected using a sterile needle.

To remove residual CsCl, the sample was dialyzed first for 2 h, and then overnight at 8°C to SM buffer (100 mM calcium chloride, 8 mM magnesium sulfate, 50 mM Tris–HCl, gelatine 0.01%), using Slide-A-Lyzer MINI Dialysis Devices (10 K MWCO; Thermo Scientific). The final purified phage preparation was subjected to titer determination via a double agar overlay plaque assay.

Purified phage particles were subjected to one-dimensional sodium dodecyl sulfate–polyacrylamide gel electrophoresis (1D SDS-PAGE) following the protocol described by Boulanger (2009). Phage protein samples were heat-denatured and subsequently separated on a 12% polyacrylamide gel (Bio-Rad). The gel was then stained using Instant Blue (Expedeon) to visualize protein bands.

A comparative analysis of the genomic regions encoding receptor-binding systems and their corresponding protein models revealed notable similarities among Dhillonvirus JLBYU60, Traversvirus tv24B, and Littlefixvirus 98PfluR60PP, as well as a correlation with Lambdavirus lambda, despite the considerable evolutionary distance between these phages (Figure 1A). Several structurally conserved domains and proteins identified in the aforementioned genomes are highlighted in Figure 1B (indicated by colored rectangles).

Lambdavirus lambda is a non-enveloped bacteriophage characterized by a long, non-contractile tail and extended lateral tail fibers. Its receptor-binding system comprises three proteins. The first, gpJ (NP_040600.1), forms the central tail spike and contains an immunoglobulin-like (Ig-like) domain at its C-terminus, which recognizes the primary receptor, the LamB protein. The second (NP_040602.1) and third (NP_040604.1) proteins constitute the tail fibers. The homo-trimer of NP_040604.1 forms the distal portion of the tail fiber and functions as a supporting receptor-binding protein. Structurally, this homo-trimer is similar to the C-terminal domain of gp37 from phage T4, which interacts with the OmpC porin. In lambda phage this region is responsible for binding both outer membrane porins C and F (Meyer et al., 2012). Between the genomic regions encoding the tail spike and the tail fibers lies an ORF encoding the outer membrane protein Lom (NP_040601.1), a β-barrel protein with an N-terminal membrane localization signal peptide. Lom is incorporated into the host membrane during infection and is expressed during lysogeny in E. coli.

Dhillonvirus JLBYU60, another non-enveloped bacteriophage with a long non-contractile tail and long lateral fibers, has been shown by Lewis et al. (2023) to use the ferrichrome outer membrane transporter (FhuA) as its main receptor. Bioinformatic analysis revealed a receptor-binding system organization comparable to that of lambda (Figure 1B). The system includes both tail fibers and a tail spike. The terminal regions of the tail fiber proteins possess carbohydrate-binding domains and intramolecular autocatalytic chaperones, suggesting an affinity for sugar motifs. Structural modelling indicated that the tail spike lacks a C-terminal Ig-like domain, instead terminating in a knife-like fold. The ORF immediately downstream encodes a small hypothetical protein (YP_010740311.1) with a β-sheet-rich structure, confirmed through structural modelling. The following ORF encodes a 23 kDa protein (YP_010740312.1), annotated by Lewis et al. as a putative membrane-associated protein. Modelling of this protein revealed a novel hand-like β-sheet fold (Figure 1B).

The next two analyzed phages, Traversvirus tv24B and Littlefixvirus 98PfluR60PP, are small non-enveloped phages with short, non-contractile tails and long lateral fibers. In Traversvirus tv24B, genome annotation showed an ORF encoding a supporting receptor-binding protein (YP_009168163.1) located ~2000 nucleotides upstream of the tail spike-encoding region. Structural analysis revealed that this protein consists of a C-terminal sugar-binding domain. The subsequent protein (YP_009168168.1) encodes a typical tail spike structurally similar to lambda’s gpJ, including N-terminal domains required for virion attachment, but lacking a C-terminal Ig-like domain—paralleling the configuration seen in Dhillonvirus JLBYU60. Structural modelling of the downstream ORF product (YP_009168169.1), annotated as an outer membrane protein, showed high structural similarity to YP_010740311.1 from Dhillonvirus. Likewise, the next downstream protein (YP_009168170.1), annotated as a membrane-associated protein, was structurally similar to YP_010740312.1 from Dhillonvirus (Figure 2). This protein, which occupies the third position in the receptor-binding region, was also annotated by the original submitters as a Lom-like outer membrane protein—an annotation confirmed by our structural bioinformatics analysis (Figure 1B).

In Littlefixvirus 98PfluR60PP, the genomic arrangement of ORFs encoding the receptor-binding system mirrors that of Dhillonvirus JLBYU60. Structural modelling showed that the tail spike protein is architecturally similar to those described above. No Ig-like receptor-binding domain was detected. Instead, downstream ORFs encode small proteins that are structurally similar to those found in Dhillonvirus and Traversvirus, although annotated here as hypothetical proteins. The final protein in the receptor-binding module (YP_010659262.1) is annotated as a minor tail protein. Structural modelling of this protein in its predicted homo-trimeric state indicates that it forms tail fibers with C-terminal carbohydrate-binding domains (Figure 1B).

Sequence and structure comparisons of these proteins are presented in Figure 2A. Despite low sequence identity, structural modelling revealed significant conservation among proteins encoded downstream of the tail spike genes.

Further comparison between the C-terminal domains of tail spikes and small Ig-like proteins encoded downstream of the spike-encoding ORFs revealed strong structural similarity to the C-terminal Ig-like domain of lambda’s gpJ (NP_040600.1), although in the analysed phages, these domains appear to be encoded as separate ORFs (Figure 2E).

Proteins encoded downstream of the Ig-like domain ORFs were generally annotated as membrane-associated proteins (represented in green in Figures 1, 2). Membrane affinity estimation (Table 1) showed that these hypothetical proteins exhibit significantly lower predicted affinity for the bacterial membrane compared to lambda’s Lom protein, suggesting that they likely do not embed into the membrane. Based on genomic localization and transcriptional direction, we hypothesize that these proteins are expressed between the tail spike and tail fiber genes and may serve structural or chaperone-like roles in tail assembly.

The obtained results suggest that the small immunoglobulin-like proteins (grey models) and hypothetical hand-like proteins (green models) may have been functionally misannotated. Based on our analyses, we hypothesized that these proteins represent previously uncharacterized structural components of the tail spike receptor-binding system. To test this hypothesis, we performed in silico simulations of tail spike assembly using heteromeric modelling with AlphaFold2-Multimer.

For Dhillonvirus JLBYU60, heteromeric modelling of the tail spike C-terminal homo-trimer in complex with the hypothetical protein YP_010740311.1 (also in a homo-trimeric form) resulted in a stable structure with a low Predicted Aligned Error (PAE) score at the interaction interface. The predicted binding affinity (ΔG) for this interaction suggested a thermodynamically favourable and potentially stable complex (Figures 3A–C).

Additional modelling of the tail spike C-terminal domain with monomeric YP_010740312.1 yielded a heteromeric structure with a higher predicted local distance difference test (pLDDT) score (Figure 3D) compared to the complex involving YP_010740311.1. Interface analysis indicated that the YP_010740312.1 monomer interacts with all three chains of the tail spike. Although the predicted binding affinity was slightly lower (Figure 3E), it was still sufficient to indicate a stable interaction. Furthermore, the PAE score for this complex was lower (Figure 3F), suggesting a more reliable structural model.

Similar modelling of the Littlefixvirus 98PfluR60PP tail spike showed comparable results to those obtained for Dhillonvirus JLBYU60 (Figure 4). The tail spike C-terminal domain was able to form stable complexes with both the homo-trimer of YP_010659264.1 and the monomer of YP_010659263.1. Notably, the interaction with YP_010659263.1 produced better pLDDT and PAE scores (Figures 4A,C,D,F) compared to YP_010659264.1. However, the predicted binding affinity was slightly higher for the homo-trimeric interaction (Figures 4B,E).

For Traversvirus tv24B, heteromeric modelling of the interaction between the tail spike C-terminal domain and YP_009168169.1 did not yield a structure resembling the immunoglobulin-like trimeric complexes observed in other phages. Nonetheless, the model predicted a localized interaction at the base of the knife-like domain, with strong binding affinity and favorable PAE scores (Figures 5A–C). Modelling of the interaction with the hand-like protein YP_009168170.1 produced better structural confidence (higher pLDDT and lower PAE scores), although the predicted binding affinity was lower (Figures 5D–F).

Altogether, modelling the spatial arrangement of the tail spike and associated hypothetical proteins did not conclusively determine the precise interaction mode. To gain further insight, we modelled interactions involving all candidate proteins for each phage. Across all analysed cases, AlphaFold2-Multimer consistently favoured tail spike termination with hand-like proteins rather than immunoglobulin-like homo-trimers. Both pLDDT and PAE scores were more favourable for the hand-like protein complexes (Figure 6), supporting our hypothesis that these structures may represent preferred or more stable configurations.

The newest bioinformatics tools for phage genome functional annotation are based on the PHROG database and use sequence or structure homology for protein function detection. Table 2 presents the comparison of functional annotation performed on the ‘newest’ phage proteins functions predictor implemented in the Sphea workflow. The ‘hand-like’ proteins from phages Dhillonvirus JLBYU60 and Traversvirus tv24B matching to PHROG 425 ‘membrane associated protein’ and the same functional annotation are assigned by pharokka, phold and phynteny. In the case of YP_010659263.1 from Littlefixvirus 98Pflur60PP the PHROG was not detected, pharokka annotated protein as ‘hypothetical’, phold and phynteny predict function as membrane associated protein.

In the case of the ‘Ig-like” small protein performed analysis matching protein YP_009168169.1from Traversvirus tv24B to PHROG 557 ‘outer membrane protein’ and functional annotated by all softwares as outer membrane protein. In the case of YP_010740311.1 from Dhillonvirus JLBYU60 and YP_010659263.1 from Littlefixvirus 98Pflur60PP matching frog were not detected, but YP_010740311.1 was functional annotated by all used software as ‘outer membrane protein’. In the case of YP_010659263.1 Pharokka did not predict protein function, but phold and phynteny annotated protein as membrane associated protein.

For Dhillonvirus JLBYU60, experimental confirmation of the main receptor was provided by Jessica M. Lewis and colleagues (Lewis et al., 2023), who identified the ferrichrome outer membrane transporter (FhuA) as the main receptor for this phage. In the present study, we performed structural modelling of interactions between candidate phage proteins and the FhuA receptor, as shown in Figures 7A–C.

Heteromeric modelling of the tail spike C-terminal domain revealed a potential interaction interface between the knife-like domain and the flexible extracellular loops of FhuA. The predicted binding affinity for this interaction was −16.9 kcal/mol. However, the modelled interface configuration suggests the parallel position of the tail spike to the bacterial cell surface, potentially limiting accessibility. The associated PAE matrix showed a high error score at the interaction site, indicating lower confidence in the predicted interface (Figure 7A).

Modelling the interaction between the homo-trimeric YP_010740311.1 and FhuA resulted in the lowest pLDDT score, the highest PAE error, and the weakest predicted binding affinity among the tested candidates (Figure 7B).

In contrast, interaction modelling between the hand-like protein YP_010740312.1 and FhuA showed a markedly higher pLDDT score and a low PAE error score across the binding interface. The predicted binding affinity indicated a strong and thermodynamically favorable interaction between YP_010740312.1 and FhuA (Figure 7C).

For Traversvirus tv24B, the receptor has been experimentally identified as the β-barrel transmembrane domain of BamA (Smith et al., 2007). Our modeling of the interaction between the tail spike C-terminal domain and BamA β-barrel domain resulted in the lowest pLDDT score, highest PAE error, and weakest predicted binding affinity of the analysed complexes (Figure 7D). Similar to Dhillonvirus, the interaction geometry positioned the receptor-binding system nearly parallel to the bacterial surface.

In contrast, modelling of the homo-trimeric YP_009168169.1 in complex with the C-terminal domain of BamA produced a low pLDDT score and high PAE error, but a high predicted binding affinity—likely due to multiple contact points across the interface (Figure 7E). The most favorable interaction was observed for YP_009168170.1: this model achieved the highest pLDDT score and the lowest PAE error at the binding interface. Although its predicted binding affinity was lower than that of the Ig-like trimer model, it remained within a thermodynamically plausible range (Figure 7F).

To evaluate the stability of the predicted protein complexes and determine whether they represent realistic assemblies rather than artifacts of computational modelling, molecular dynamics (MD) simulations were performed. The simulations were conducted in an aqueous environment at 300 K over a 1 μs timescale.

Root-mean-square deviation (RMSD) analysis (Figure 8A) indicated that the complex formed between the tail spike C-terminal domain and the hand-like protein remained more stable throughout the simulation compared to the complex with the Ig-like protein. The latter underwent significant structural rearrangements within the first 100 ns of the simulation, with deviations exceeding 10 Å, before reaching a more stabilized state. In contrast, the RMSD of the hand-like protein complex remained around 5 Å throughout the full 1 μs simulation.

Root-mean-square fluctuation (RMSF) analysis further supported these findings. The tail spike C-terminal domain complexed with the Ig-like protein exhibited greater overall flexibility, particularly in the peripheral regions of the Ig-like domains, which showed fluctuations exceeding 7 Å (Figure 8B). Conversely, the complex with the hand-like protein displayed significantly lower RMSF values, indicating greater structural rigidity. The two regions of highest flexibility in this complex (RMSF > 5 Å) corresponded to external loop regions of the hand-like protein (Figure 8C).

A video of the MD simulation for both structural layouts is available in Supplementary material 2.

Protein modelling results suggest that the tail spikes of the analysed phages may be terminated by an additional protein exhibiting high affinity for receptors on the bacterial cell surface. However, these in silico predictions require experimental validation.

The phage sample after ultracentrifugation is presented in Figure 9A, the lower grey bands were collected and analysed. To assess the correlation between theoretically identified structural proteins and those isolated from virions, SDS-PAGE analysis was performed for Littlefixvirus 98PfluR60PP (Figure 9B). Figure 9C presents the digital reconstruction of the analysed virions. This phage was originally described by Kazimierczak and colleagues as an effective and safe agent against pathogenic, antibiotic-resistant Pseudomonas spp. in the study (Kazimierczak et al., 2019). The transmission electron microscopy picture of this bacteriophage is available in this thesis.

The SDS-PAGE results were compared with the theoretical molecular masses of structural proteins predicted through functional annotation, allowing for the tentative identification of the observed bands (Table 3). A prominent band at approximately 21 kDa corresponded well with the hypothetical protein YP_010659263.1—identified in this study as a candidate hand-like protein potentially acting as a missing receptor-binding element.

Additionally, a diffused band was observed around 13 kDa, which may correspond to the hypothetical Ig-like protein or alternatively to a minor capsid or scaffolding protein. The intensity of this band suggests a relatively high abundance of the protein in the virion, supporting the possibility that it plays a structural role, potentially as a minor capsid component or internal scaffold.