In silico characterisation of lysis and host-recognition modules in Staphylococcus aureus bacteriophage genomesAbstract
BACKGROUND Methicillin-resistant Staphylococcus aureus (MRSA) is a leading human pathogen, and declining antibiotic efficacy has renewed interest in lytic bacteriophages, whose efficacy phenotypic assays predict poorly.
OBJECTIVES To characterise the lysis and host-recognition modules of lytic S. aureus bacteriophages, and to determine how far the detectability of each module depends on the annotation source rather than on the genomes.
METHODS Twenty-two S. aureus phage genomes from NCBI GenBank were processed through an identical detection funnel under three annotation sources: the deposited record, Pharokka and Phold. Modules were validated by InterPro, disrupted endolysins resolved by tBLASTn, and each domain aligned separately to measure identity.
FINDINGS The LysK-type endolysin (CHAP-Amidase-SH3b) was conserved and detected identically under all three sources. Tail fibre/receptor-binding protein (RBP) detection was strongly source-dependent: 14 of 22 genomes from the deposited annotation, 20 of 22 after re-annotation, and 22 of 22 combined. Resolved by domain, conservation was scale-dependent: the amidase domain was invariant within species while SH3b varied by 5.7%.
MAIN CONCLUSIONS Lysis-module conservation is independent of the annotation source; the apparent variability of host recognition is largely an artefact of incomplete deposited annotation. Module-level surveys require re-annotation before comparative claims.
Key words:
bacteriophages;
Staphylococcus aureus
; methicillin-resistant Staphylococcus aureus; genome; viral; phylogeny; host specificity
Staphylococcus aureus, and in particular methicillin-resistant S. aureus (MRSA), is among the most clinically important human bacterial pathogens, responsible for infections that range from skin and soft-tissue infections to life-threatening bacteraemia, endocarditis, and medical device-associated disease, and it is a high-priority target for the development of new antibacterial therapeutics.(1) Antimicrobial resistance (AMR) driven by MRSA has escalated into a major global health threat, intensifying the need for precision non-antibiotic therapeutics.(1,2) The structural resilience of MRSA, mediated by the mecA gene product PBP2a and by OatA-dependent peptidoglycan O-acetylation, undermines the efficacy of both conventional and last-line antibiotics.(3,4) Lytic bacteriophages offer a rational alternative, delivering precise bactericidal activity through their lysis cassette — the holin and endolysin proteins that degrade peptidoglycan from within the host cell — without disrupting the commensal microbiota.(5,6) Infection begins one step earlier, at the host-recognition cassette. The receptor-binding proteins (RBPs) of a tailed phage are the tail fibres, tailspikes and baseplate-associated proteins that make first contact with the bacterial surface; in Staphylococcus these recognise wall teichoic acid and its glycosylation pattern. RBPs are modular: an N-terminal region anchors the protein to the virion while a C-terminal region carries the receptor-binding specificity, so that the host range of a phage can be redirected by exchanging or mutating the C-terminal region alone, without altering the lytic enzymes.(7,8) Lysis and host recognition are therefore complementary determinants of a phage’s therapeutic behaviour, and the present study characterises both.
Despite the conceptual elegance of phage lysis, phenotypic evaluation in vitro often gives unpredictable results, because several biological barriers intervene.(9,10) A paradoxical observation illustrates this limitation: phage vB_SauP-436A markedly reduced MRSA biofilm biomass, whereas phage vB_SauM-515A1 increased biofilm biomass, even though both phages exhibited high planktonic lytic efficiency in vitro.(11,12) Lytic efficacy alone is therefore an insufficient predictor of therapeutic success, particularly in chronic biofilm-associated infection.(13)
Complete genome sequences of S. aureus bacteriophages are now abundantly available in public repositories such as NCBI GenBank, and most tailed phages have recently been reclassified within the class Caudoviricetes, including the families Herelleviridae and Rountreeviridae, superseding the earlier morphology-based grouping.(14) The product names in these deposited records, however, are often inconsistent with the underlying protein domain architecture. Re-examining the lysis and host-recognition modules against curated domain databases therefore offers an efficient and reproducible way to establish their functional content without isolating new phages.(15,16)
Although many S. aureus phage genomes are publicly available, the lysis and host-recognition modules of therapeutically relevant lytic phages have not been compared systematically. This study therefore set out to characterise the lysis and host-recognition modules of lytic S. aureus bacteriophages and to map candidate genomic features that may underlie inter-phage differences in therapeutic response. Because the dataset is limited to publicly deposited sequences, the study is descriptive and comparative. It is intended to provide a genomic foundation for selecting natural phage candidates and to inform later experimental work, not to dictate clinical formulation.
MATERIALS AND METHODS
Genome retrieval and dataset curation - Genome sequences were retrieved from NCBI GenBank between 18 and 30 May 2026 using the search query ‘Staphylococcus phage complete genome’, supplemented by targeted retrieval of the specific strains examined in a previous study(12) and downloaded as GenBank flat files. The initial set was curated against predefined criteria. The inclusion criteria were (i) records annotated as complete genome sequences, (ii) a lytic lifestyle inferred from genome annotation and the absence of lysogeny-associated genes (integrase, excisionase, CI repressor), and (iii) S. aureus as the annotated host. The exclusion criteria were (i) duplicate or near-identical genome sequences, (ii) the presence of lysogeny-associated genes, and (iii) jumbo genomes exceeding 200 kb. During curation, Staphylococcus phage vB_SauP_EBHT (NC_055906) was excluded because its intergenomic similarity to Staphylococcus phage Portland (MT926124) is 100.0%, and Staphylococcus phage MarsHill (MW248466) was excluded for exceeding 200 kb. The retained dataset was deliberately not dereplicated at the species level: because the object of study is the consistency of annotation rather than the extent of genomic diversity, genome pairs above the species demarcation threshold are informative rather than redundant, and their retention is what makes the annotation discordance reported below detectable. Staphylococcus phage Sb1M_6168 (MN336262) and Sb1M_9832 (MN336263) were retained on that basis as representatives of the Sb-1 lineage; the full similarity matrix is given in Supplementary data (Table I). One record, Staphylococcus phage vB_SauM-515A1 (MN047438), carried a Draft/Partial status flag and was retained as a single, explicitly noted exception to the complete-genome criterion because all modules analysed here are fully annotated and the phage is of direct comparative relevance to the previous study. Host identity of all 22 retained records was manually confirmed against the GenBank source entries prior to analysis. The final dataset comprised 22 genomes.
Intergenomic similarity - Nucleotide-level relatedness among the retained genomes, together with the two genomes excluded during curation, was quantified with VIRIDIC (Virus Intergenomic Distance Calculator)(17) using the default BLASTN parameter set (-word_size 7 -reward 2 -penalty -3 -gapopen 5 -gapextend 2) and the demarcation thresholds recommended by the ICTV Bacterial and Archaeal Viruses Subcommittee, namely 95% intergenomic similarity for species and 70% for genus. Aligned genome fractions and genome length ratios were recorded alongside each similarity value, since a high similarity supported by a low aligned fraction does not indicate genome-wide relatedness.
Genomic feature extraction - Genome metadata were extracted using custom scripts written in Python 3.12.10 with Biopython 1.87(18) and pandas 3.0.3. Using the Entrez and SeqIO modules, the scripts parsed each record to obtain genome size, GC content, the number of coding sequences (CDS), and the number of tRNA genes. Taxonomic ranks (class, family, subfamily) were resolved from the taxonomy field of each record by ICTV rank suffix (-viricetes, -viridae, -virinae); subfamilies placed by NCBI in no family (the Azeredovirinae phages EW and SA13) were recorded as family Unassigned. CDS counts are reported as deposited in GenBank; de novo re-annotation by independent gene callers may yield slightly different totals (for example, 236 versus 238 predicted features for MN047438). Transfer RNA counts in Table I are likewise those of the deposited record; they are reported separately from the de novo predictions described below, and the two are not interchangeable. Re-annotation and annotation-source comparison. All extracted values were verified manually against the source records and are summarised in Table I.
Functional annotation of lysis and host-recognition modules - Endolysin-candidate proteins were extracted from each genome and submitted to InterPro Version 108.0(19) for domain architecture analysis against the Pfam, CDD, SUPERFAMILY, and Gene3D databases; this identified the CHAP (cysteine, histidine-dependent amidohydrolase/peptidase), Amidase, SH3b (bacterial Src-homology 3b), LysM, and NlpC/P60 domains reported in Table II. Holin and tail fibre/receptor-binding protein (RBP) genes were detected from CDS product annotations using curated keyword sets, with each positive call retained together with the matched annotation for verification; the structural tail-tube protein, present in essentially all tailed phages, was deliberately excluded from the RBP keyword set. Because this procedure depends on the wording of the deposited product name, whereas endolysins were additionally validated at the domain level, the two module classes are not supported by evidence of equal strength. This asymmetry was addressed by re-annotation (see below), and RBP candidates were in addition submitted to InterPro under the same protocol used for endolysins, so that the evidence type supporting every call in Table II is stated explicitly. Where the automated candidate did not correspond to a recognised peptidoglycan hydrolase (Maine) or no intact endolysin-sized open reading frame was present (Sb1M_6168 and Sb1M_9832), the endolysin was treated as not confidently assigned by automated extraction; for the Sb-1 lineage the catalytic fragment was recovered by tBLASTn using the Sb1_8383 endolysin as query. In two further genomes, vB_SauM-515A1 (MN047438) and vB_SauM-fRuSau02 (MF398190), the endolysin (lysK) is annotated as one moiety because the gene is interrupted by a self-splicing intron; the intact enzyme is reconstituted at the RNA level, as demonstrated transcriptomically for vB_SauM-515A1,(20) and the corresponding Table II entries are annotated accordingly. Targeted tBLASTn searches used the Sb1_8383 endolysin (MN336261, 495 aa) as query against each genome with BLAST+ 2.16.0, -outfmt 6 and -max_target_seqs 5, all other parameters at their defaults (e-value 10, word size 3, BLOSUM62, gap open 11, gap extend 1, SEG filtering on); a hit was accepted as the genuine endolysin at 90% identity or above. The outcome of each search is recorded in the host profile as a curated case with its identity value and literature source, so that re-running the pipeline reproduces the resolved assignment without repeating the search. Throughout, an endolysin is described as intact when its catalytic and cell-wall-binding domain architecture is complete in the mature protein, whether or not that protein is encoded by a single uninterrupted open reading frame; the term refers to domain architecture and not to sequence identity, which is reported separately in Supplementary data (Table II).
Domain architecture of lysis modules and presence of host-recognition modules among the 22 bacteriophages
Domain-partitioned pairwise identity - To quantify what “conserved” refers to, the unique endolysin sequences were partitioned at their InterPro-resolved domain boundaries and each domain set was aligned separately, so that variation in one domain could not be averaged away by another. Boundaries were taken from the InterPro output already generated for Table II rather than re-predicted; where several member databases annotated the same domain, the widest span was used and its source recorded (CHAP, PS50911; amidase, G3DSA:3.40.80.10; SH3b, G3DSA:2.30.30.40, the Pfam amidase signature PF01510 being returned for too few sequences to define the set). Each partition, and the full-length sequences, were aligned on the MAFFT online service with the L-INS-i strategy(20) and pairwise percent identity was computed over columns in which both sequences carry a residue, columns with a gap in either sequence being excluded from numerator and denominator alike. Because the three domains were not resolved in the same set of proteins, identity was compared across domains only within the subset carrying all three.
Phylogenetic analysis - Phylogenetic relationships were inferred from the large terminase subunit (TerL), a conserved marker suitable for resolving relationships across phage families.(21,22) TerL sequences were extracted from each genome; an exact-match check was required to capture the abbreviated ‘Ter’ product name used in the Kayvirus group. Where a genome carried more than one terminase-annotated coding sequence, the longest was taken, since terminase genes are frequently split across several open reading frames and returning the first match recovers a fragment while the full-length product lies further along the genome. This affected Twort (NC_007021), whose terminase is interrupted by HNH endonuclease insertions and annotated as three separate coding sequences; the longest fragment (286 aa against a 605-aa reference) was used, InterPro returned no domain for it, and the position of Twort in Figure therefore rests on partial data. Its separation from the remaining genomes is independently supported by intergenomic similarity (34.5-36.8%, below the ICTV genus demarcation threshold). Staphylococcus phage Portland (MT926124) and vB_SauP-436A1 (MN150710) lacked an annotated TerL — consistent with their atypical small genome sizes (< 20 kb) — and were excluded, leaving 20 sequences. Sequences were aligned on the MAFFT online service using the L-INS-i strategy with all other settings at their defaults.(20) A maximum-likelihood tree was constructed in MEGA 12.1.2(23) under the LG+G+I substitution model, with rate variation modelled by a discrete Gamma distribution and a proportion of invariant sites; node support was assessed with 1000 bootstrap replicates, gaps were treated by partial deletion with an 80% site-coverage cutoff, and the tree was rooted on the Azeredovirinae lineage, represented here by Staphylococcus phages EW (NC_007056) and SA13 (NC_021863), which are recovered as a pair at 100% bootstrap and are the only genomes in the dataset placed by NCBI in no family. The tree was annotated in iTOL with five tip-level layers drawn from Tables I, II and III: family, genome size, endolysin domain architecture, lysis-gene disruption (HNH insertion and intron split), and the receptor-binding protein call before and after re-annotation.
Detection of holin and receptor-binding protein (RBP) across three annotation sources in 22 lytic Staphylococcus aureus bacteriophage genomes
Maximum-likelihood phylogenetic tree of 20 Staphylococcus aureus bacteriophages based on large terminase subunit (TerL) sequences, inferred in MEGA 12.1.2 under the LG+G+I model with 1000 bootstrap replicates, rooted on the Azeredovirinae lineage (EW, NC_007056 and SA13, NC_021863) and annotated in iTOL. Tip labels give the phage name; accessions are listed in Table I. Tip annotations, from left to right: the receptor-binding protein call, where the grey square is the call from the deposited annotation and the orange square the call after re-annotation; family; genome size in kb; endolysin domain architecture; and lysis-gene disruption, where a filled circle marks an HNH endonuclease insertion and a filled star an intron-split lysK. The eight genomes carrying an empty grey square beside a filled orange one are those whose receptor-binding protein was recovered only after re-annotation, and they fall within the group whose members are at or above the species demarcation threshold to one another. Bootstrap values are shown at the nodes: support is 100% for the 17-phage group and for the EW-SA13 pair, and 27-64% for nodes within the 17-phage group, where internal branch lengths approach zero. An asterisk marks Twort, represented by the longest of the three coding sequences into which its terminase is split (286 aa). Portland (MT926124) and vB_SauP-436A1 (MN150710) carry no terminase-annotated coding sequence and are not included. Scale bar, substitutions per site.
Re-annotation and annotation-source comparison - To test whether module detection reflects genomic content or annotation completeness, all 22 genomes were re-annotated independently of their deposited records and re-processed through the identical detection funnel, so that the annotation layer was the only variable. Structural gene calling and functional assignment were performed with Pharokka v1.10.1;(24) the resulting GenBank records were passed to Phold v1.3.0,(25) which assigns function by ProstT5 protein-language-model embeddings followed by Foldseek structural search against the Phold database, executed on an NVIDIA RTX 4060 GPU at half precision (PyTorch 2.8.0+cu128); residual proteins still annotated as hypothetical after Phold were passed to Phynteny_transformer v0.1.3, which assigns PHROG (Prokaryotic Virus Remote Homologous Groups) functional categories from genomic context using ESM-2 embeddings. Phynteny returns a PHROG category rather than a product name and was therefore not used to make module calls; its output is reported separately in Supplementary data (Table III). Three detection arms were compared: the deposited GenBank annotation, the Pharokka re-annotation, and the Phold re-annotation. All three used the same keyword sets, the same InterPro protocol and the same curation funnel; only the source of the CDS product names differed. Results are given in Table III.
Codon usage and phage-encoded tRNAs - Relative synonymous codon usage (RSCU)(26) was computed for the host by pooling the coding sequences of five curated S. aureus reference genomes (NCTC 8325, NC_007795.1; Mu50, NC_002758.2; MRSA252, NC_002952.2; USA300_FPR3757, NC_007793.1; and CP023390.1) rather than a single strain, and, separately, for each phage genome from its complete CDS complement, using a single shared implementation so that the definition could not differ between host and phage. Similarity between each phage and the host was quantified over the 59 codons belonging to degenerate synonymous families (methionine and tryptophan excluded, being invariant at RSCU = 1) as the Pearson correlation and the Spearman rank correlation of the two RSCU vectors, together with GC content at the third codon position (GC3) and the mean absolute RSCU difference. Transfer RNA genes were predicted de novo with tRNAscan-SE 2.0(27) rather than read from the deposited annotation, and each predicted anticodon was mapped to the codons it decodes and to the host RSCU of those codons; a codon was classified as rare in the host at RSCU < 0.5.
Ethics - This study is entirely computational and analysed only publicly available genome sequences; no ethics approval was required and no biological material was accessed.
RESULTS
General genomic characteristics of lytic S. aureus bacteriophages - The 22 complete bacteriophage genomes retrieved from NCBI GenBank varied substantially (Table I). Genome size ranged more than eight-fold, from 17,471 bp in Staphylococcus phage Portland (MT926124) to 148,564 bp in Staphylococcus phage 676Z (NC_047724). GC content varied within a narrower range, from 29.34% in vB_SauP-436A1 (MN150710) to 35.99% in phage EW (NC_007056), and the number of predicted CDS ranged from 20 in phage Portland to 236 in both vB_SauM-fRuSau02 (MF398190) and vB_SauM-515A1 (MN047438). On the basis of genome size and CDS count, the phages fell into three size classes: a micro class (17-18 kb, 20-21 CDS) comprising Portland and vB_SauP-436A1, both assigned to the family Rountreeviridae; a medium class (42-45 kb, 62-77 CDS) comprising phages EW and SA13; and a macro class (127-148 kb, 164-236 CDS) encompassing the majority of the dataset, including the Kayvirus group and the Sb-1 lineage. GC content was tightly constrained: 20 of the 22 phages fell below 31%, slightly below the 32-33% genomic GC of S. aureus. Two phages, EW (35.99%) and SA13 (35.42%), were clear outliers with GC content above the host average. Transfer RNA genes were annotated in only five of the deposited records (four tRNAs each in vB_SauM-fRuSau02, A5W, vB_SauM-515A1, Sb-1 and Maine), whereas de novo scanning recovered tRNA genes in 18 of the 22 genomes: three tRNAs (Asp-GTC, Ile2-CAT, Phe-GAA) in 17 macro-class genomes and a single Ile2-CAT in Twort. The four genomes returning no tRNA de novo were the two micro-class genomes, Portland and vB_SauP-436A1, and the two medium-class genomes, EW and SA13. Per-genome counts from the de novo scan are provisional: the calibration control recovered three of its four expected tRNAs, so the absolute number per genome may be underestimated by one, although the presence/absence pattern across genomes is unaffected.
Lysis modules and the source-dependence of host-recognition detection - Functional annotation of the 22 genomes revealed high conservation of the intracellular lysis components, and an apparent variability of the extracellular host-recognition components that the following section shows to be a property of the deposited annotation rather than of the genomes (Table II). Among phages with a canonical free endolysin, protein length clustered tightly around 480-495 amino acids, and fourteen of these displayed the classic LysK-type architecture, defined by tandem N-terminal catalytic domains (CHAP and Amidase) and a C-terminal SH3b cell wall-binding domain (CBD); proteins falling outside this range corresponded to specific structural or genomic exceptions rather than to true size variation of the endolysin. In the deposited annotation, holin was detected in 21 of the 22 phages and tail fibre or RBP genes in 14 of the 22. As shown in the following section, both figures are properties of the deposited annotation rather than of the genomes, and neither should be read as a measure of biological variability. The conserved LysK-type architecture extended even to the phylogenetically divergent phage SA13 (481 amino acids), which retained the full CHAP-Amidase-SH3b organisation despite its long branch and elevated GC content; the principal exception was phage EW, the most divergent lineage, whose 482-amino-acid endolysin diverges from the canonical LysK architecture.
Two mechanistically distinct departures from the canonical architecture were observed, both of which confound keyword-based annotation. In the Sb-1 lineage (Sb1M_6168, Sb1M_9832) the endolysin open reading frame is interrupted by an HNH endonuclease insertion, leaving no intact endolysin-sized frame in the deposited inventory; in vB_SauM-515A1 and vB_SauM-fRuSau02 the lysK gene is interrupted by a self-splicing group I intron carrying the ksaI homing endonuclease, and the intact CHAP-Amidase-SH3b enzyme is reconstituted at the RNA level.(28,29) In three further genomes the keyword extractor returned a protein that was not the free endolysin - a non-lytic YbiA-like N-glycosidase in Maine, and virion-associated peptidoglycan hydrolases in JD007 and Twort - and a targeted tBLASTn search using the Sb1_8383 endolysin as query was required to resolve each case, recovering a canonical 495-amino-acid LysK at 99% identity in Maine and JD007 but only a partial, divergent region in Twort.(30) The individual cases are described in Supplementary data.
Conservation is reported here as pairwise amino-acid identity within a defined domain, not as the mere presence of a domain, since a domain can be present in every protein while its sequence varies widely. Partitioning the endolysins on that basis shows that “conserved” applies at two different scales and with opposite rank order. Six of the unique sequences carried all three domains and could be compared directly. Among the five drawn from macro-class genomes, which are themselves at or above the species demarcation threshold [Supplementary data (Table I)], the amidase domain was invariant at 100% identity across all ten pairs, CHAP ranged from 99.2 to 100%, and SH3b from 94.3 to 100%, making SH3b the only domain to fall below 95% at this scale. Between SA13 and those five sequences the order reverses: CHAP is the most divergent domain at 22.9% identity in all five pairs, SH3b intermediate at 47.6-48.6%, and the amidase domain the least divergent at 54.2-54.5%, against 42.7-43.0% for the full-length proteins. The distribution is strictly bimodal, with ten pairs above 94% and five below 60% and none between, so a single summary statistic across all pairs would describe neither population [Supplementary data (Table II)].
Detection of lysis and host-recognition modules depends on the annotation source - Re-processing the same 22 genomes through the identical funnel under three annotation sources changed the detection of both host-recognition and holin genes, while leaving endolysin detection unaffected (Table III). Holin rose from 21 of 22 in the deposited annotation to 22 of 22 under both Pharokka and Phold. Tail fibre/RBP genes rose from 14 of 22 to 20 of 22 under both re-annotations, and the total number of RBP loci rose from 26 to 57 (Pharokka) and 58 (Phold). Taking the union of the three sources, holin and RBP were each detected in 22 of 22 genomes; no genome was negative in all three. The change is concentrated in the Sb-1 lineage and the Kayvirus group. Sb1_8383, Sb1M_6168 and Sb1M_9832 each moved from no detected RBP locus to three, while Sb-1 itself remained at three, removing a discordance that had appeared within a single lineage. Five further macro-class genomes (676Z, Fi200W, MSA6, P4W, qdsa002) moved from zero to three loci.
Two genomes moved in the opposite direction - vB_SauP-436A1 and Portland, the two micro-class Rountreeviridae genomes, were RBP-positive in the deposited annotation but negative under both re-annotations. Inspection of the underlying calls indicates that the product names assigned by Pharokka and Phold to the short tail fibres of these genomes fall outside our keyword set; this is a limitation of the detection layer applied here, not a biological difference, and it is reported as such. Both genomes are positive in the union column. Endolysin detection was invariant across the three sources: an endolysin representative was selected in the same 18 genomes in every arm, collapsing to the same 11 unique sequences and the domain architectures reported in Table II were unchanged. InterPro domain support was recovered for the RBP candidates of 11 of 22 genomes in each arm, although not for the same 11: eight genomes were supported in all three arms, and the union across arms was 14 of 22. In every arm the domain-supported figure is lower than the keyword-based call, consistent with the sparse representation of receptor-binding protein folds in profile-HMM databases and with structural homology being the more appropriate comparator for this module class. Passing the Phold records through Phynteny and re-processing the resulting 22 genomes through the same funnel changed nothing: holin was detected in 22 of 22 genomes, tail fibre/RBP in 20 of 22, and the number of RBP loci was 58, matching the Phold arm genome for genome in all three measures. Phynteny is therefore reported as the residual annotation it resolves [Supplementary data (Table III)] rather than as a fourth detection arm.
Intergenomic similarity provides an annotation-independent test of this interpretation - Across the 22 retained genomes, 67 of the 231 genome pairs reach or exceed the ICTV species demarcation threshold of 95%, and three pairs reach 100.0%: A5W with Staph1N, Sb1M_6168 with Sb1M_9832, and the excluded Portland with vB_SauP_EBHT. Of those 67 species-level pairs, 38 are discordant in the deposited RBP call, one genome being scored positive and the other negative, the most extreme case being vB_SauM-fRuSau02 and MSA6 at 99.8% intergenomic similarity. Genome pairs that are 99-100% identical at the nucleotide level cannot differ in the presence or absence of an entire gene family, so these discordances are necessarily artefacts of the deposited annotation rather than genuine differences in gene content. Because this test uses no re-annotation of any kind, it corroborates Table III independently.
At the genus demarcation threshold of 70%, the dataset resolves into six clusters: a single macro-class cluster of 17 genomes; a Rountreeviridae cluster containing vB_SauP-436A1, Portland and the excluded vB_SauP_EBHT; and Twort, EW and SA13 each alone. Twort shares only 34.5-36.8% intergenomic similarity with the macro-class cluster despite its basal position within the TerL core clade, and EW and SA13 share 27.4% with each other and no more than 0.5% with any other genome, consistent with their long branches and elevated GC content [Supplementary data (Table I)].
Phylogenetic relationships based on the large terminase subunit - The maximum-likelihood tree inferred from TerL sequences, rooted on phage EW, separated the 20 analysed phages into a supported core group, a basal lineage and a divergent pair (Figure). Seventeen phages with large genomes (127-148 kb) consistent with a myovirus-type morphology formed a group supported at 100% bootstrap, within which internal branch lengths were near zero and internal node support ranged from 27 to 64%. Staphylococcus phage Twort (NC_007021) fell outside that group as its immediate sister lineage, consistent with its recognised status as a prototype myovirus and with its intergenomic similarity to the group (34.5-36.8%) lying below the genus demarcation threshold. Phages EW (NC_007056) and SA13 (NC_021863) formed the divergent pair, joined at 100% bootstrap, with terminal branches of 0.797 and 1.545 substitutions per site respectively and a shared internal branch of 1.961. Of the 190 pairwise TerL distances, 92 were exactly zero. The tree therefore resolves lineage separation but not branching order within the 17-phage group, and Figure should be read accordingly.
Codon usage similarity to the host is stratified by size class and coincides with tRNA carriage - Codon usage of every phage was compared with that of S. aureus [Supplementary data (Table IV)]. Similarity was high overall but not uniform, the Pearson correlation of RSCU vectors ranging from 0.844 to 0.972 (Spearman 0.797-0.956). The dataset separated completely into two groups with no intermediate value. The four genomes carrying no de novo-detectable tRNA — vB_SauP-436A1, Portland, EW and SA13 — showed the closest correspondence to host codon usage (r = 0.927-0.972). The 18 genomes carrying tRNAs showed the lowest (r = 0.844-0.881), the two ranges being separated by a gap between 0.881 and 0.927.
This ordering is not reducible to base composition. Host GC3 was 22.73%; the two micro-class genomes lay closest to it (21.28% and 21.31%), but EW and SA13, whose GC3 of 30.62% and 29.85% departs from the host by more than seven percentage points, nonetheless resembled the host more closely in RSCU than the macro-class genomes, whose GC3 of 18.70-20.37% lies within four points of the host value. Codon preference within synonymous families therefore carries information that GC3 alone does not. The largest single contributions to the macro-class divergence were an excess of AGA against a deficit of CGT among the arginine codons, and an excess of CCT against a deficit of CCA among the proline codons.
The tRNA complement of the macro-class genomes was uniform: Asp-GTC, Ile2-CAT and Phe-GAA, decoding five codons, of which one (GAC) is rare in the host at RSCU < 0.5. In every case the remaining members of the same synonymous family were also decoded, so the complement does not constitute selective supplementation of rare host codons; it is equally consistent with generic translational support and with the proposed anti-defence role, and the present sequence-based analysis does not distinguish between them.
DISCUSSION
Genomic characterisation stratifies these phages in a way that phenotypic data alone cannot. The earlier experimental evaluation of lytic activity(12) did not treat genomic capacity as a variable influencing therapeutic outcome. Yet the more than eleven-fold difference in CDS count between Portland (20 CDS) and vB_SauM-fRuSau02 (236 CDS) implies large differences in coding capacity for lysis strategy, accessory functions and engineering potential, even though the accessory-gene repertoires behind those differences were not annotated here. This stratification is consistent with established morphological descriptions, in which larger contractile-tailed (myovirus-type) phages correspond to the macro class and smaller phages to the micro class; macro-class phages possess sufficient coding capacity to accommodate additional engineering modules, an advantage absent in micro-class phages whose genomes approach the functional minimum.(31)
The conserved LysK-type architecture acts cooperatively to cleave the pentaglycine cross-bridges characteristic of the S. aureus peptidoglycan,(32) and its conservation even in the divergent phage SA13 underscores the strong selective constraint on the lytic enzyme. Against this conserved background, the two independent mechanisms of gene interruption identified here — HNH endonuclease insertion in the Sb-1 lineage and self-splicing intron interruption of lysK in the Kayvirus group — are notable chiefly because they mislead annotation that relies on product names alone: in both cases the intact catalytic architecture is preserved, and the apparent anomaly is an artefact of the deposited annotation rather than a functional deficiency. Reliable module assignment therefore requires domain-level validation and homology search rather than on textual product names, a methodological point with direct consequences for the curation of phage genomes used in therapeutic candidate selection. Measured directly, that conservation must be qualified in one respect. Among endolysins from genomes of the same species, the catalytic domains are effectively fixed, the amidase domain being identical in every pair and CHAP differing by less than 1%, whereas the cell-wall-binding SH3b domain varies by up to 5.7%. Recognition-level variation is therefore not confined to the tail module: it is also present, at a smaller amplitude, in the binding domain of an otherwise invariant lytic enzyme. This is consistent with reports that SH3b-type binding domains of staphylococcal endolysins differ in strain-binding profile, and with the modular engineering of phage scaffolds in which binding specificity is exchanged independently of catalytic function,(8) and it means that the conservation of the lysis module should be read as conservation of catalytic architecture rather than of the entire molecule. Between lineages the ranking reverses, CHAP being the most divergent domain, so no single domain can be described as the conserved or the variable one without stating the comparison scale.
The contrast between a conserved lysis module and an apparently variable host-recognition module does not survive re-annotation. Under domain- and structure-level annotation, host-recognition genes are detectable in every genome examined, and the variability visible in the deposited records is a property of those records. What the comparison does establish is that the two module classes differ in how reliably they are annotated: endolysin assignments were stable across all three sources, whereas holin and RBP assignments were not. This interpretation bears directly on the paradoxical biofilm behaviour of vB_SauM-515A1 reported previously:(12) because its endolysin is intron-split rather than truncated, and the intact LysK retains the full catalytic architecture, its atypical biofilm response is unlikely to stem from endolysin deficiency; host recognition, translational compatibility, matrix-acting depolymerases and anti-defence systems remain equally plausible and as-yet-untested contributors. The two phages whose contrasting biofilm behaviour motivated this survey differ on a further axis that the present analysis makes measurable. vB_SauP-436A1 has the codon usage most similar to that of S. aureus in the dataset (r = 0.972) and encodes no tRNA, whereas vB_SauM-515A1 has the least similar (r = 0.844) and carries the three-tRNA macro-class complement; they also differ more than eight-fold in genome size and coding capacity. These are properties of the genome sequence, not measurements of biofilm activity, and no causal inference is drawn from them; they are noted because they identify translational compatibility, alongside host recognition, as a candidate variable that the original phenotypic comparison did not control and that would be testable experimentally. The paradox that motivated this survey is not confined to Staphylococcus. Phage predation has been reported to increase rather than reduce biofilm biomass across several bacterial genera, through mechanisms as different as lysis-driven release of extracellular DNA in Shewanella oneidensis,(33) structural incorporation of filamentous phage into the Pseudomonas aeruginosa matrix,(34) and coevolutionary selection for mucoid phenotypes in Pseudomonas fluorescens.(35) A broader survey of biofilm formation under phage predation reaches the same conclusion.(36) That the same paradoxical direction arises from unrelated mechanisms suggests it reflects a general property of phage-biofilm interaction rather than a peculiarity of any one host, and reinforces the central limitation of a genomic survey such as this one: sequence capacity does not predict the direction of the biofilm response. Whether the reconstituted enzyme is kinetically equivalent to a canonical LysK is not established. Splicing restores the CHAP-Amidase-SH3b architecture at the RNA level and the mature protein is reported as full length,(28) so no structural basis for reduced activity is apparent; but neither specific activity nor turnover has been measured for the spliced product, and any residual difference would also depend on splicing efficiency in vivo, which is likewise unmeasured. The consequence for the present interpretation is the same either way. If the reconstituted enzyme is equivalent, the atypical biofilm response of vB_SauM-515A1 must be attributed to factors outside the lysis module; if it is not, the difference is a kinetic property of the enzyme rather than the absence or truncation of a module, and would still be invisible to a genomic survey. Distinguishing the two requires purified-enzyme kinetics on the spliced product, which lies outside the scope of a sequence-based analysis.
The phylogenetic analysis reinforces this conclusion. The near-zero internal branch lengths and consequent low internal support within the 17-phage group reflect the high conservation of TerL among these phages rather than a methodological limitation, since a single conserved marker such as TerL resolves inter-family relationships well but discriminates poorly among closely related intra-family members.(21) Ninety-two of the 190 pairwise TerL distances were exactly zero, which is what that conservation amounts to in practice. The long branches and elevated GC content of EW and SA13 are mutually consistent and support independent evolutionary trajectories for these two phages. Importantly, the placement of vB_SauM-515A1 within the conserved core clade indicates that its atypical biofilm behaviour cannot be attributed to deep phylogenetic divergence. Nucleotide-level comparison refines this picture: Twort, although basal to the TerL 17-phage group, falls below the ICTV genus threshold against every other macro-class genome (34.5-36.8%), so its basal position on a single-marker tree reflects the conservation of the marker rather than membership of the same genus.
Among the accessory features catalogued here, phage-encoded tRNAs are of particular interest, and their distribution is informative only when the deposited annotation is set aside. De novo prediction recovered tRNA genes in 18 of 22 genomes rather than the five in which they had been annotated, and the four genomes without them are precisely the four whose codon usage most closely matches that of the host. Carriage of tRNAs and divergence from host codon usage therefore coincide completely in this dataset, with no genome falling between the two groups. The direction of that relationship cannot be established here — a phage whose codon usage has drifted from its host may be under selection to carry its own tRNAs, or carriage may relax the constraint that would otherwise keep its codon usage host-like — and the dataset contains no independent measure of translational efficiency with which to test either. Phage-encoded tRNAs have been proposed both as translational supplements that support phage protein synthesis when the host tRNA pool is suboptimal and, more recently, as anti-defence elements that counteract bacterial immune systems;(37,38,39) the composition observed here does not discriminate between these roles, since the single host-rare codon decoded (GAC) is accompanied by decoding of the rest of its synonymous family. Both roles remain hypotheses requiring experimental confirmation in S. aureus phages and were not tested in the present, sequence-based analysis. More generally, the lysis and host-recognition cassettes are not the only modules that could contribute to differences in therapeutic response. The translational module examined here, comprising phage-encoded tRNAs and the codon usage of the phage genome, varies measurably across the dataset and separates it completely into two groups, and other modules not examined here, such as anti-defence systems and depolymerases acting on the biofilm matrix, are equally plausible contributors. The present survey characterises two cassettes because they are the two for which domain-level evidence can be obtained from sequence alone; it does not imply that they are the only relevant ones. This bears on what the survey can and cannot be used for. It does not license a selection strategy that simply seeks phages whose lysis and host-recognition modules most closely resemble the canonical ones, and the present results argue against such a strategy: the canonical architecture is so nearly universal in this dataset that it does not discriminate between candidates, and the two phages with the most divergent lysis modules are not the two with the divergent biofilm phenotypes. What a module-level survey does provide is the opposite service: it identifies which genomes carry a disrupted or non-canonical module and therefore require experimental verification before use, and it provides the domain-level inventory against which a phenotypic result can later be interpreted. Candidate selection remains a phenotypic decision informed by genomics, not replaced by it.
Several limitations should be acknowledged. The analysis is annotation-dependent: feature counts and module assignments derive from deposited GenBank records and computational domain inference, which approximate but do not replace experimental validation of protein function. The magnitude of that dependence is quantified in Table III: re-annotation changed the RBP call in eight of 22 genomes and the holin call in one, without any change to the underlying sequences. Detection also depends on the keyword vocabulary applied downstream of annotation, as the two micro-class genomes demonstrate; module surveys that report presence or absence without stating both the annotation source and the detection vocabulary are not comparable between studies. The study is descriptive and comparative; no statistical hypothesis testing was applied, and the modest size of the dataset limits generalisation. The phylogeny rests on a single conserved marker and therefore does not capture the mosaic, recombination-rich nature of phage genomes. Accordingly, the genomic features identified here are best regarded as hypothesis-generating signals that require subsequent experimental confirmation, in line with the role of this study as preliminary groundwork for planned experimental investigation. The codon usage comparison is descriptive and correlational: RSCU similarity and tRNA carriage covary completely in this dataset, but with 22 genomes falling into two size-stratified groups the association cannot be separated from lineage, and no causal or directional claim is made. The domain-partitioned identity analysis rests on the six unique endolysins in which all three domains were resolved, five of which come from genomes above the species demarcation threshold; the within-species figures therefore describe one lineage rather than the dataset as a whole. The phylogeny rests on a single conserved marker in genomes that are largely conspecific, so it supports statements about lineage separation but not about branching order within the 17-phage group; and the Twort tip is a partial sequence, as stated in Methods. No enzymatic measurement was made: where the manuscript reports an intact domain architecture it reports genomic capacity, not catalytic equivalence to a canonical enzyme.
In conclusion - This comparative genomic survey of lytic S. aureus bacteriophages shows that the lysis module is conserved at the level of domain architecture, and that this conservation is recovered identically from deposited and re-annotated records. Host-recognition genes appeared variable in the deposited annotation, being detectable in 14 of 22 genomes, but were recovered in 20 of 22 after re-annotation and in all 22 when the annotation sources were combined; their apparent variability therefore reflects annotation completeness rather than genomic content. Domain-level analysis distinguished true endolysins from virion-associated peptidoglycan hydrolases and revealed two independent mechanisms by which lysis-module genes can be interrupted — HNH endonuclease insertion and intron splitting — both of which can mislead annotation that relies on product names alone. The same effect extends beyond the lysis and host-recognition modules: de novo prediction recovered tRNA genes in 18 of 22 genomes against five in the deposited records, and the genomes lacking them proved to be those whose codon usage most closely matches the host. Module-level comparisons of phage genomes therefore require the annotation source to be stated, and re-derived wherever comparative claims are made. These findings provide a genomic foundation for the rational selection of phage candidates and a basis for experimentally investigating planktonic and biofilm-associated responses in subsequent work, in which S. aureus, as a major human pathogen, remains the clinically relevant target.
SUPPLEMENTARY MATERIALS
ACKNOWLEDGEMENTS
To the Department of Biology, Faculty of Science and Technology, Universitas Terbuka.
How to cite:
DATA AVAILABILITY
All genome sequences analysed in this study are publicly available from the NCBI GenBank database under the accession numbers listed in Table I. The GenBank flat files (.gb) used as input, together with the custom Python analysis scripts that reproduce all reported results, are openly available in a public repository (https://github.com/HivanA98/In-silico-characterisation-of-lysis-and-host-recognition-modules-in-S.-aureus-bacteriophage-genomes), archived at a fixed version under release tag v3.2.3-MIOC-revision and deposited with the persistent identifier https://doi.org/10.5281/zenodo.21988477. No new nucleotide or protein sequences were generated in this work; therefore, no additional public-bank accession numbers were created. The re-annotated GenBank records produced by Pharokka, Phold and Phynteny, the three-arm detection outputs underlying Table III, and the codon usage matrices underlying Supplementary data (Table IV) are deposited in the same repository under a tagged release.
REFERENCES
-
1 WHO - World Health Organization. Global antimicrobial resistance and use surveillance system (GLASS) report 2022. Geneva: WHO; 2022. https://iris.who.int/handle/10665/364996
» https://iris.who.int/handle/10665/364996 -
2 Nimmana BK, Nguyen AD. Antibiotic resistance. In: StatPearls. Treasure Island: StatPearls Publishing; 2026. Available from: https://www.ncbi.nlm.nih.gov/books/NBK513277/
» https://www.ncbi.nlm.nih.gov/books/NBK513277/ -
3 Jones CS, Sychantha D, Howell PL, Clarke AJ. Structural basis for the O-acetyltransferase function of the extracytoplasmic domain of OatA from Staphylococcus aureus J Biol Chem. 2020; 295(24): 8204-13. https://doi.org/10.1074/jbc.RA120.013108.
» https://doi.org/10.1074/jbc.RA120.013108 -
4 Sari MAR, Qurrohman MT, Dewangga VS. Detection of mecA gene as a marker for Staphylococcus aureus types of methicillin-resistant Staphylococcus aureus (MRSA) using PCR technique. J Biol Tropis. 2024; 24(4): 204-11. https://doi.org/10.29303/jbt.v24i4.7539
» https://doi.org/10.29303/jbt.v24i4.7539 -
5 Brüser T, Mehner-Breitfeld D. Occurrence and potential mechanism of holin-mediated non-lytic protein translocation in bacteria. Microb Cell. 2022; 9(10): 159-73. https://doi.org/10.15698/mic2022.10.785.
» https://doi.org/10.15698/mic2022.10.785 -
6 Liu K, Wang C, Zhou X, Guo X, Yang Y, Liu W, et al. Bacteriophage therapy for drug-resistant Staphylococcus aureus infections. Front Cell Infect Microbiol. 2024; 14: 1336821. https://doi.org/10.3389/fcimb.2024.1336821.
» https://doi.org/10.3389/fcimb.2024.1336821 -
7 Yehl K, Lemire S, Yang AC, Ando H, Mimee M, Torres MDT, et al. Engineering phage host-range and suppressing bacterial resistance through phage tail fibre mutagenesis. Cell. 2019; 179(2): 459-469.e9. https://doi.org/10.1016/j.cell.2019.09.015.
» https://doi.org/10.1016/j.cell.2019.09.015 -
8 Ando H, Lemire S, Pires DP, Lu TK. Engineering modular viral scaffolds for targeted bacterial population editing. Cell Syst. 2015; 1(3): 187-96. https://doi.org/10.1016/j.cels.2015.08.013.
» https://doi.org/10.1016/j.cels.2015.08.013 -
9 Duarte AC, Fernández L, De Maesschalck V, Gutiérrez D, Campelo AB, Briers Y, et al. Synergistic action of phage phiIPLA-RODI and lytic protein CHAPSH3b: a combination strategy to target Staphylococcus aureus biofilms. NPJ Biofilms Microbiomes. 2021; 7(1): 39. https://doi.org/10.1038/s41522-021-00208-5.
» https://doi.org/10.1038/s41522-021-00208-5 -
10 Anyaegbunam NJ, Anekpo CC, Anyaegbunam ZKG, Doowuese Y, Chinaka CB, Odo OJ, et al. The resurgence of phage-based therapy in the era of increasing antibiotic resistance: from research progress to challenges and prospects. Microbiol Res. 2022; 264: 127155. https://doi.org/10.1016/j.micres.2022.127155.
» https://doi.org/10.1016/j.micres.2022.127155 -
11 Abdraimova NK, Shitikov EA, Malakhova MV, Gorodnichev RB, Kornienko MA. Effects of lytic bacteriophages of the families Herelleviridae and Rountreeviridae on the Staphylococcus aureus biofilms. Bull Russ State Med Univ. 2024; (6): 98-104. https://doi.org/10.24075/brsmu.2024.061. Also available from: https://vestnik.rsmu.press/archive/2024/6
» https://doi.org/10.24075/brsmu.2024.061» https://vestnik.rsmu.press/archive/2024/6 -
12 Hasugian IA, Zikriyani H. Literature review on the therapeutic potential of bacteriophages against resistant Staphylococcus aureus J Biol Tropis. 2025; 25(4): 4869-79. https://doi.org/10.29303/jbt.v25i4.9938.
» https://doi.org/10.29303/jbt.v25i4.9938 -
13 Molendijk MM, Boekema BKHL, Lattwein KR, Vlig M, Bode LGM, Koopmans MPG, et al. Bacteriophage therapy reduces Staphylococcus aureus in a porcine and human ex vivo burn wound infection model. Antimicrob Agents Chemother. 2024; 68(9): e00650-24. https://doi.org/10.1128/aac.00650-24.
» https://doi.org/10.1128/aac.00650-24 -
14 Barylski J, Kropinski AM, Alikhan NF, Adriaenssens EM. ICTV virus taxonomy profile: Herelleviridae. J Gen Virol. 2020; 101(4): 362-3. https://doi.org/10.1099/jgv.0.001392.
» https://doi.org/10.1099/jgv.0.001392 -
15 Hatfull GF, Hendrix RW. Bacteriophages and their genomes. Curr Opin Virol. 2011; 1(4): 298-303. https://doi.org/10.1016/j.coviro.2011.06.009.
» https://doi.org/10.1016/j.coviro.2011.06.009 -
16 Moller AG, Petit RA, Read TD. Species-scale genomic analysis of Staphylococcus aureus genes influencing phage host range and their relationships to virulence and antibiotic resistance genes. mSystems. 2022; 7(2): e01083-21. https://doi.org/10.1128/msystems.01083-21.
» https://doi.org/10.1128/msystems.01083-21 -
17 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.
» https://doi.org/10.3390/v12111268 -
18 Cock PJA, Antao T, Chang JT, Chapman BA, Cox CJ, Dalke A, et al. Biopython: freely available Python tools for computational molecular biology and bioinformatics. Bioinformatics. 2009; 25(11): 1422-3. https://doi.org/10.1093/bioinformatics/btp163.
» https://doi.org/10.1093/bioinformatics/btp163 -
19 Blum M, Andreeva A, Florentino LC, Chuguransky SR, Grego T, Hobbs E, et al. InterPro: the protein sequence classification resource in 2025. Nucleic Acids Res. 2025; 53(D1): D444-56. https://doi.org/10.1093/nar/gkae1082.
» https://doi.org/10.1093/nar/gkae1082 -
20 Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019; 20(4): 1160-6. https://doi.org/10.1093/bib/bbx108.
» https://doi.org/10.1093/bib/bbx108 -
21 Meier-Kolthoff JP, Göker M. VICTOR: genome-based phylogeny and classification of prokaryotic viruses. Bioinformatics. 2017; 33(21): 3396-404. https://doi.org/10.1093/bioinformatics/btx440.
» https://doi.org/10.1093/bioinformatics/btx440 -
22 Dion MB, Oechslin F, Moineau S. Phage diversity, genomics and phylogeny. Nat Rev Microbiol. 2020; 18(3): 125-38. https://doi.org/10.1038/s41579-019-0311-5.
» https://doi.org/10.1038/s41579-019-0311-5 -
23 Stecher G, Suleski M, Tao Q, Tamura K, Kumar S. MEGA 12.1: cross-platform release for macOS and Linux operating systems. J Mol Evol. 2026; 94(1): 14-8. https://doi.org/10.1007/s00239-025-10287-z.
» https://doi.org/10.1007/s00239-025-10287-z -
24 Bouras G, Nepal R, Houtak G, Psaltis AJ, Wormald PJ, Vreugde S. Pharokka: a fast scalable bacteriophage annotation tool. Bioinformatics. 2023; 39(1): btac776. https://doi.org/10.1093/bioinformatics/btac776.
» https://doi.org/10.1093/bioinformatics/btac776 -
25 Bouras G, Grigson SR, Mirdita M, Heinzinger M, Papudeshi B, Mallawaarachchi V, et al. Protein structure informed bacteriophage genome annotation with Phold. Nucleic Acids Res. 2026; 54(1): gkaf1448. https://doi.org/10.1093/nar/gkaf1448.
» https://doi.org/10.1093/nar/gkaf1448 -
26 Sharp PM, Tuohy TMF, Mosurski KR. Codon usage in yeast: cluster analysis clearly differentiates highly and lowly expressed genes. Nucleic Acids Res. 1986; 14(13): 5125-43. https://doi.org/10.1093/nar/14.13.5125.
» https://doi.org/10.1093/nar/14.13.5125 -
27 Chan PP, Lin BY, Mak AJ, Lowe TM. tRNAscan-SE 2.0: improved detection and functional classification of transfer RNA genes. Nucleic Acids Res. 2021; 49(16): 9077-96. https://doi.org/10.1093/nar/gkab688.
» https://doi.org/10.1093/nar/gkab688 -
28 Kornienko M, Fisunov G, Bespiatykh D, Kuptsov N, Gorodnichev R, Klimina K, et al. Transcriptional landscape of Staphylococcus aureus Kayvirus bacteriophage vB_SauM-515A1. Viruses. 2020; 12(11): 1320. https://doi.org/10.3390/v12111320.
» https://doi.org/10.3390/v12111320 -
29 Kornienko M, Bespiatykh D, Gorodnichev R, Abdraimova N, Shitikov E. Transcriptional landscapes of Herelleviridae bacteriophages and Staphylococcus aureus during phage infection: an overview. Viruses. 2023; 15(7): 1427. https://doi.org/10.3390/v15071427.
» https://doi.org/10.3390/v15071427 -
30 Latka A, Maciejewska B, Majkowska-Skrobek G, Briers Y, Drulis-Kawa Z. Bacteriophage-encoded virion-associated enzymes to overcome the carbohydrate barriers during the infection process. Appl Microbiol Biotechnol. 2017; 101(8): 3103-19. https://doi.org/10.1007/s00253-017-8224-6.
» https://doi.org/10.1007/s00253-017-8224-6 -
31 Pires DP, Meneses L, Brandão AC, Azeredo J. An overview of the current state of phage therapy for the treatment of biofilm-related infections. Curr Opin Virol. 2022; 53: 101209. https://doi.org/10.1016/j.coviro.2022.101209.
» https://doi.org/10.1016/j.coviro.2022.101209 -
32 Gutiérrez D, Rodríguez-Rubio L, Ruas-Madiedo P, Fernández L, Campelo AB, Briers Y, et al. Design and selection of engineered lytic proteins with Staphylococcus aureus decolonizing activity. Front Microbiol. 2021; 12: 723834. https://doi.org/10.3389/fmicb.2021.723834.
» https://doi.org/10.3389/fmicb.2021.723834 -
33 Godeke J, Paul K, Lassak J, Thormann KM. Phage-induced lysis enhances biofilm formation in Shewanella oneidensis MR-1. ISME J. 2011; 5(4): 613-26. https://doi.org/10.1038/ismej.2010.153.
» https://doi.org/10.1038/ismej.2010.153 -
34 Secor PR, Sweere JM, Michaels LA, Malkovskiy AV, Lazzareschi D, Katznelson E, et al. Filamentous bacteriophage promote biofilm assembly and function. Cell Host Microbe. 2015; 18(5): 549-59. https://doi.org/10.1016/j.chom.2015.10.013.
» https://doi.org/10.1016/j.chom.2015.10.013 -
35 Scanlan PD, Buckling A. Co-evolution with lytic phage selects for the mucoid phenotype of Pseudomonas fluorescens SBW25. ISME J. 2012; 6(6): 1148-58. https://doi.org/10.1038/ismej.2011.174.
» https://doi.org/10.1038/ismej.2011.174 -
36 Hosseinidoust Z, Tufenkji N, van de Ven TGM. Formation of biofilms under phage predation: considerations concerning a biofilm increase. Biofouling. 2013; 29(4): 457-68. https://doi.org/10.1080/08927014.2013.779370.
» https://doi.org/10.1080/08927014.2013.779370 -
37 Azam AH, Kondo K, Chihara K, Nakamura T, Ojima S, Nie W, et al. Evasion of antiviral bacterial immunity by phage tRNAs. Nat Commun. 2024; 15(1): 9586. https://doi.org/10.1038/s41467-024-53789-y.
» https://doi.org/10.1038/s41467-024-53789-y -
38 Lomeli-Ortega CO, Balcázar JL. Why tRNA acquisition could be relevant to bacteriophages? Microb Biotechnol. 2024; 17(4): e14464. https://doi.org/10.1111/1751-7915.14464.
» https://doi.org/10.1111/1751-7915.14464 -
39 van den Berg DF, Brouns SJJ. Phage tRNAs: decoding the enigma. Trends Microbiol. 2025; 33(10): 1121-31. https://doi.org/10.1016/j.tim.2025.05.009.
» https://doi.org/10.1016/j.tim.2025.05.009
Correspondence
+ Corresponding author: hivanarmadi@gmail.comConflict of interest
Handling editor:
FIRST REVIEW ROUND - REVIEWERS’ COMMENTS
REVIEWER #1
General Comments
Overall, the manuscript addresses an interesting topic by exploring the lysis and host-recognition modules in S. aureus phages to guide non-traditional therapeutic choices. While the study is primarily descriptive, the authors adequately acknowledge its limitations. My main concern relates to what is defined as “conserved” throughout the text, as this concept remains unclear and requires further specification.
Specific Points
The English language could be improved.
Lines 12-18 (Page 3): Has this type of paradox only been observed in Staphylococcus, or does it also occur in other species? A brief mention of this broader context would enrich the discussion.
Could only the lysis and host-recognition modules reflect differences in therapeutic response, or might other modules also play a role? Consequently, does this imply that the strategy for identifying clinically relevant phages relies on finding those with lysis and host-recognition modules similar to canonical ones? Please clarify this rationale.
The following reference is inaccessible: Abdraimova NK 2024 (doi: 10.24075/brsmu.2024.061). Please verify and provide a working link or correct the citation.
Is the codon usage of these phages similar to that of S. aureus? Phage replication is expected to be more efficient in a host with a similar codon usage profile. Discussing this could add depth to the host-phage dynamics.
In the Introduction, the lysis cassettes (holin and endolysin) are mentioned, but what constitutes the host-recognition cassette? Furthermore, since the study verifies the functional domains of these modules, it would be beneficial to briefly introduce these domains in the Introduction as well.
In Table II, it would be helpful to visually or categorically differentiate between the lysis and the host-recognition modules.
Lines 40-42 (Page 6) and Table II: What exactly does the “conservation” of these components refer to? Is it presence/absence, amino acid sequence length, or the presence of specific domains? It would be highly relevant to report the sequence identity of these sequences. Even if functional domains are conserved, sequence variations in other regions might indicate distinct selective pressures.
Lines 31-32 (Page 9): Again, does “intact” refer to the protein’s length or to 100% sequence identity? Please specify.
Lines 54-55 (Page 6): Please define the RBP genes. Additionally, could you briefly explain how the tail fiber determines the host range?
In the results, lysins are well described, detailing various cases where LysK is non-canonical. However, to improve the manuscript’s readability and dynamism, I would suggest moving a portion of these descriptive details to a Supplementary Material.
I expected a more specific comparative analysis of the phages vB_SauP-436A and vB_SauM-515A1, which were cited in the Introduction as inducing distinct biofilm biomasses. Elaborating on this comparison would strengthen the study’s narrative connection to the introduction.
Line 46 (Page 9): References 24 and 26 do not seem appropriate to support this statement.
Lines 47-48 (Page 9): Would an intron-split endolysin, once processed, have the same catalytic efficiency as a canonical endolysin? Could there be structural or kinetic differences influencing the biofilm response? If the efficiency is indeed the same, then the biofilm’s behavior must be attributed to other factors.
REVIEWER #2
Well-curated dataset, explicit criteria, clear writing. The most valuable contribution is re-examining deposited annotations against curated domain databases instead of trusting product names. That result currently sits as an intermediate step, while the result given prominence is unsupported.
Major
1. The central comparison is not supported and the paper needs reframing around what it actually shows. No host-range data are analysed, yet the Objectives and Conclusion assert that host recognition governs host range and is the rational target for engineering. The asymmetry beneath this claim is the core problem: endolysins were characterised with InterPro, tBLASTn and manual curation, while holin and RBP were scored by keyword matching alone. Domain-level characterisation is being compared against text matching and the difference reported as biology. The 14/22 figure cannot mean what it is taken to mean — a tailed phage cannot lack a receptor-binding apparatus, and Table II shows the artefact directly, with Sb-1 scored positive while Sb1_8383, Sb1M_6168 and Sb1M_9832 are negative despite being treated as one lineage. Holin fails the same way. Both figures are lower bounds on annotation completeness, not measures of variability.
2. Re-annotate with Pharokka, Phold and Phynteny. The thesis is that deposited annotations are unreliable, yet the entire analysis rests on them. This is what makes comment 1 answerable rather than merely conceded.
3. “Highly conserved” and “markedly variable” carry no measurement. Supply pairwise identity across the canonical endolysins resolved by domain — if SH3b varies while the catalytic domains do not, recognition is not confined to the tail module and the central dichotomy collapses. The same imprecision governs a curation decision: one genome was excluded as “essentially identical” to another while two genomes of identical length were retained. State the threshold.
4. Reproducibility. Missing: keyword sets, tBLASTn parameters, MAFFT version, curation funnel, and a DOI with fixed commit. Deposit the TerL alignment.
Minor
1. Table II is titled as validated by InterPro, but InterPro was applied only to endolysins; the RBP column is a binary keyword call.
2. Figure 1 does not work as a visual. Beyond resolution, the tree is currently decorative: the manuscript narrates lineages, size classes, architecture exceptions and the two disruption mechanisms in text, none of which is visible on the tree. Consider rebuilding it in iTOL or FigTree and annotating the tips with those layers, such as family, genome size, endolysin architecture, HNH insertion and intron split, and the RBP call, so the figure carries the comparison rather than restating it.
3. The Ethics section is disproportionate for a computational study.
4. The AI-assistance declaration is misplaced. Its content is appropriate, but Acknowledgements is for human and institutional contributions that do not meet authorship criteria, and a tool cannot be acknowledged in that sense. Move it to a dedicated Declaration of AI Use section, or to Conflict of Interest if the journal has no such heading.
AUTHORS’ RESPONSE TO THE REVIEWERS
Dear Dr Vicente,
We thank you and both reviewers for a review that improved this manuscript substantially. Reviewer 2 identified a claim that our data did not support, and correcting it required re-analysing the entire dataset under three annotation sources rather than one. That work changed the paper’s central result, and we are grateful for it. Reviewer 1’s insistence on defining “conserved” led us to measure what we had previously only asserted, with a result that partly contradicts our original framing.
Every comment is addressed below in the order received. Changes in the manuscript are highlighted. All code, genome records, alignments, intermediate results and the host profile are deposited at https://doi.org/10.5281/zenodo.21988477 under release tag v3.2.3-MIOC-revision.
Yours sincerely,
Ivan Armadi Hasugian, on behalf of both authors
Summary of principal changes
Reframing — No claim about host range remains. The 14/22 figure is reported as a lower bound on annotation completeness.
New Table III — Detection of holin and RBP across three annotation sources, all 22 genomes.
New Supplementary Table S1 — Residual annotation resolved by Phynteny.
New Supplementary Table S2 — Phage-versus-host codon usage with de novo tRNA counts.
New Supplementary Table S3 — VIRIDIC intergenomic similarity, 24 genomes.
New Supplementary Table S4 — Endolysin per-domain pairwise identity.
New Supplementary Text S1 — Case-by-case description of non-canonical endolysin assignments.
Figure 1 — Rebuilt in iTOL with five annotation layers; tree corrected after a fix to the terminase extractor; supplied as TIFF at 300 dpi.
Table II — New RBP InterPro domain column, corrected caption, two new footnotes, host-recognition columns shaded.
New Methods subsections — Intergenomic similarity; Domain-partitioned pairwise identity; Re-annotation and annotation-source comparison; Codon usage and phage-encoded tRNAs.
New references — 31-39, covering RSCU, tRNAscan-SE, Pharokka, Phold, VIRIDIC and four cross-genus biofilm studies.
Data availability — Zenodo DOI 10.5281/zenodo.21988477, release tag v3.2.3-MIOC-revision.
Reviewer 2
We have placed Reviewer 2 first because the reframing requested in Major 1 determined the structure of the whole revision.
Major 1 — The central comparison is not supported and the paper needs reframing
Comment. No host-range data are analysed, yet the Objectives and Conclusion assert that host recognition governs host range. Endolysins were characterised with InterPro, tBLASTn and manual curation, while holin and RBP were scored by keyword matching alone. The 14/22 figure cannot mean what it is taken to mean.
Response. We accept this in full. The asymmetry was real and the claim was unsupported. We re-processed all 22 genomes through the identical S1-S4 detection funnel under three annotation sources - the deposited GenBank records, Pharokka v1.10.1 and Phold v1.3.0 (ProstT5 embeddings with Foldseek structural search) - so that the annotation layer was the only variable. Holin detection rose from 21/22 to 22/22 and tail fibre/RBP from 14/22 to 20/22 under both re-annotations, with RBP loci rising from 26 to 57 and 58. The four genomes the reviewer named behave exactly as predicted: Sb1_8383, Sb1M_6168 and Sb1M_9832 each move from zero to three RBP loci while Sb-1 remains at three, removing the within-lineage discordance. Taking the union of the three sources, holin and RBP are each detected in 22/22 genomes.
Response. The analysis requested in Major 3 supplied an independent confirmation that does not depend on re-annotation at all. Of the 231 genome pairs among the 22 retained genomes, 67 reach or exceed the ICTV species demarcation threshold of 95%, and 38 of those 67 are discordant in the deposited RBP call - the most extreme at 99.8% intergenomic similarity. Genome pairs that are 99-100% identical cannot differ in the presence of an entire gene family, so those discordances are necessarily annotation artefacts.
Response. We also report two genomes moving in the opposite direction. vB_SauP-436A1 and Portland are RBP-positive in the deposited annotation and negative after re-annotation, because our own keyword vocabulary does not capture the product names assigned to the short tail fibres of these micro-class genomes. We state this as a limitation of the detection layer rather than a biological finding.
Action. The Objectives, Results, Discussion and Conclusion were rewritten; no claim about host range remains anywhere in the manuscript. The figure of 14/22 is now reported as a lower bound on the completeness of deposited annotation. Table III (new) gives the three-arm comparison for all 22 genomes. A new Methods subsection, “Re-annotation and annotation-source comparison”, describes the procedure, and a new Results subsection reports it. The same effect proved to extend to a third feature class and is now also reported: de novo prediction recovers tRNA genes in 18 of 22 genomes against five in the deposited records.
Major 2 — Re-annotate with Pharokka, Phold and Phynteny
Comment. The thesis is that deposited annotations are unreliable, yet the entire analysis rests on them.
Response. All three tools were run. Pharokka and Phold are reported as full detection arms in Table III. Phynteny assigns a PHROG functional category to proteins that remain of unknown function rather than a product name, so it cannot generate module calls on the same basis as the other two; presenting it as a fourth column would imply an independent confirmation that it does not provide. To settle the point empirically rather than by argument, we nevertheless split the Phynteny output into one file per genome and re-processed all 22 through the identical funnel: the results match the Phold arm genome for genome on holin presence, on RBP presence and on RBP locus count. Phynteny assigned a category to all 2,599 proteins that remained of unknown function after Phold, but at low confidence (median 0.23; 2.5% at or above 0.8), and one protein in the entire dataset was assigned to the tail category.
Action. Pharokka and Phold results are in Table III and in the new Results subsection. The Phynteny result is reported in Supplementary Table S1, with the reasoning above stated in Methods. All four arms are deposited in the code archive.
Major 3 — “Highly conserved” and “markedly variable” carry no measurement
Comment. Supply pairwise identity across the canonical endolysins resolved by domain - if SH3b varies while the catalytic domains do not, recognition is not confined to the tail module. The same imprecision governs a curation decision: one genome was excluded as “essentially identical” to another while two genomes of identical length were retained. State the threshold.
Response. Both halves are now measured, and the first partly contradicts our original framing, which we report rather than soften. The unique endolysins were partitioned at their InterPro-resolved domain boundaries and each domain set aligned separately. Among the six sequences carrying all three domains, conservation proved to be scale-dependent with opposite rank order at the two scales. Within species the amidase domain is invariant at 100% across all ten pairs, CHAP ranges 99.2-100%, and SH3b ranges 94.3-100% - so SH3b is the only domain to fall below 95%, exactly as the reviewer anticipated. Between SA13 and those five sequences the order reverses: CHAP is the most divergent at 22.9%, SH3b intermediate at 47.6-48.6%, and amidase the least divergent at 54.2-54.5%. Recognition-level variation is therefore not confined to the tail module: it is present, at smaller amplitude, in the binding domain of an otherwise invariant lytic enzyme.
Response. For the curation threshold, we quantified intergenomic similarity with VIRIDIC using the default BLASTN parameter set and the ICTV demarcation thresholds. vB_SauP_EBHT is 100.0% similar to Portland, which is why it was excluded. However, applying the same criterion consistently shows that our original rule was not applied consistently: A5W/Staph1N and Sb1M_6168/Sb1M_9832 also reach 100.0% and were retained. Rather than defend an inconsistent rule we have changed the stated rule. The dataset is deliberately not dereplicated at the species level, because the object of study is the consistency of annotation rather than the extent of genomic diversity - and it is precisely the retention of near-identical genomes that makes the annotation discordance reported under Major 1 detectable.
Action. A new Methods subsection, “Domain-partitioned pairwise identity”, describes the procedure; the results are in the Results and in Supplementary Table S4. A new Methods subsection, “Intergenomic similarity”, states the VIRIDIC parameters and thresholds; the curation paragraph was rewritten with the numerical threshold, and the full matrix is Supplementary Table S3. “Conserved” is now defined operationally in the Results as pairwise amino-acid identity within a defined domain, not as the presence of a domain.
Major 4 — Reproducibility
Comment. Missing: keyword sets, tBLASTn parameters, MAFFT version, curation funnel, and a DOI with fixed commit. Deposit the TerL alignment.
Response. All five are now supplied. The keyword sets are in the deposited host profile: holin (5 terms), RBP (12), lysis (16, a deliberately broad collection net) and an 18-term precision denylist of broad-keyword false positives confirmed non-endolysin by InterPro. The tBLASTn command, its parameters and the 90% acceptance threshold are stated in Methods; we also note there that the outcome of each search is recorded in the profile as a curated case, so that re-running the pipeline reproduces the resolved assignment without repeating the search. MAFFT was the online service, run 16 August 2026 with the L-INS-i strategy and all other settings at default; the full MEGA 12.1.2 parameter set is given in Methods. The curation funnel is documented in the archive.
Action. The archive is deposited at https://doi.org/10.5281/zenodo.21988477 under release tag v3.2.3-MIOC-revision, which fixes the commit. It contains the code, the host profile, the deposited and re-annotated genome records, the InterPro output, the VIRIDIC matrix, the per-stage results for every arm, and the alignments - including the TerL alignment and the four endolysin domain alignments. The Data Availability statement cites the DOI and the tag.
Minor 1 — Table II is titled as validated by InterPro
Comment. InterPro was applied only to endolysins; the RBP column is a binary keyword call.
Response. Correct. In the revision InterPro was additionally applied to the RBP candidates in all three annotation arms. Domain support was recovered for 11 of 22 genomes in each arm, although not for the same 11: eight were supported in all three and the union across arms was 14 of 22. This is lower than the keyword-based call in every arm, and we report it as such rather than as an improvement: receptor-binding proteins are sparsely represented in profile-HMM databases, which is why structural homology as implemented in Phold, rather than domain assignment, is the appropriate symmetric comparator for this module class.
Action. The caption of Table II was corrected to state which validation applies to which column. A new column, RBP InterPro domain, was added, and two footnotes now state the evidence type behind each call. The three host-recognition columns are shaded to separate them visually from the lysis columns, which also addresses a request from Reviewer 1. We also note, for reproducibility, that the archive contains two versions of the endolysin extraction script: S4_endolysin_extractor_for_interpro.py is the standalone script used for the original submission and S4_endolysin_extractor.py the profile-driven version used throughout the revision. Both are deposited so either version of Table II can be reproduced; the revision used the latter throughout.
Minor 2 — Figure 1 does not work as a visual
Comment. The tree is currently decorative. Consider rebuilding it in iTOL or FigTree and annotating the tips with those layers, so the figure carries the comparison rather than restating it.
Response. Figure 1 was rebuilt and now carries five tip-level layers drawn from Tables I, II and III: family, genome size, endolysin domain architecture, lysis-gene disruption (HNH insertion and intron split), and the receptor-binding protein call before and after re-annotation. The last layer makes the central result of the revision visible directly: eight genomes show an empty square for the deposited call beside a filled square for the re-annotated call, and all eight lie inside the group whose members are at or above the species demarcation threshold to one another.
Response. Rebuilding the figure also required correcting the tree itself. The terminase extractor now returns the longest terminase-annotated coding sequence rather than the first, which changes the sequence used for Twort, whose terminase is split across three coding sequences by HNH insertions. Branch lengths, node support and the rooting were updated accordingly. We also state in Methods that the Twort tip is a 286-amino-acid fragment for which InterPro returns no domain, and that its separation from the remaining genomes is independently supported by intergenomic similarity (34.5-36.8%, below the ICTV genus threshold).
Action. Figure 1 was rebuilt in iTOL, is supplied as TIFF at 300 dpi, and the legend now states the annotation layers, the bootstrap support range for nodes inside the 17-phage group (27-64%), and the partial status of the Twort sequence. The iTOL annotation files are deposited with the tree.
Minor 3 — The Ethics section is disproportionate
Comment. The Ethics section is disproportionate for a computational study.
Response. Agreed. We have reduced it to a single sentence rather than removing it altogether, so that the absence of an ethics requirement is still stated explicitly for a study that uses only public sequence data.
Action. The Ethics paragraph in Materials and Methods was reduced from a full paragraph to one sentence.
Minor 4 — The AI-assistance declaration is misplaced
Comment. Acknowledgements is for human and institutional contributions; a tool cannot be acknowledged in that sense. Move it to a dedicated Declaration of AI Use section.
Response. Agreed.
Action. The declaration was moved to a dedicated Declaration of AI Use section, and the Acknowledgements now contain only institutional acknowledgement.
Reviewer 1
General comment — the definition of “conserved”
Comment. My main concern relates to what is defined as “conserved” throughout the text, as this concept remains unclear and requires further specification.
Response. This was the most consequential comment we received from either reviewer, and answering it changed a claim rather than only a wording. Conservation is now defined operationally and measured. The full response is given under Reviewer 2, Major 3; in summary, conservation is reported as pairwise amino-acid identity within a defined domain, not as the presence of a domain, and it proved to be scale-dependent with opposite rank order within and between lineages.
Action. A definition sentence was added to the Results; the measurement is reported there and in Supplementary Table S4; a new Methods subsection describes the procedure.
- The English language could be improved
Response. The manuscript has been revised throughout for language. Overlong sentences in the Introduction, Results and Discussion were divided, redundant constructions removed, and spelling and terminology made consistent (British spelling throughout; “tail fibre” and “behaviour” standardised).
Action. Language revisions are distributed throughout and are highlighted with all other changes.
- Has this type of paradox only been observed in Staphylococcus?
Comment. Has this type of paradox only been observed in Staphylococcus, or does it also occur in other species? A brief mention of this broader context would enrich the discussion.
Response. It is not confined to Staphylococcus, and the mechanisms differ between the reported cases, which is itself informative. Phage predation has been reported to increase rather than reduce biofilm biomass through lysis-driven release of extracellular DNA in Shewanella oneidensis, structural incorporation of filamentous phage into the Pseudomonas aeruginosa matrix, and coevolutionary selection for mucoid phenotypes in Pseudomonas fluorescens. That the same paradoxical direction arises from unrelated mechanisms suggests a general property of phage-biofilm interaction rather than a peculiarity of one host.
Action. Two sentences and four references were added to the Discussion, and used to reinforce the limitation that sequence capacity does not predict the direction of the biofilm response. We kept this brief, as the reviewer asked for a mention of broader context rather than a review.
- Might other modules also play a role, and does this imply a selection strategy?
Comment. Could only the lysis and host-recognition modules reflect differences in therapeutic response, or might other modules also play a role? Consequently, does this imply that the strategy for identifying clinically relevant phages relies on finding those with lysis and host-recognition modules similar to canonical ones? Please clarify this rationale.
Response. Other modules certainly may. The translational module examined here - phage-encoded tRNAs and the codon usage of the phage genome - varies measurably across the dataset and separates it completely into two groups, and modules we did not examine, such as anti-defence systems and depolymerases acting on the biofilm matrix, are equally plausible contributors. We characterised these two cassettes because they are the two for which domain-level evidence can be obtained from sequence alone, not because they are the only relevant ones.
Response. On the second half of the question, our results argue against such a selection strategy. The canonical architecture is so nearly universal in this dataset that it does not discriminate between candidates, and the two phages with the most divergent lysis modules are not the two with the divergent biofilm phenotypes. What a module-level survey provides is the opposite service: it identifies which genomes carry a disrupted or non-canonical module and therefore require experimental verification before use, and it supplies the domain-level inventory against which a later phenotypic result can be interpreted. Candidate selection remains a phenotypic decision informed by genomics, not replaced by it.
Action. Both points were added to the Discussion.
- The Abdraimova NK 2024 reference is inaccessible
Comment. Please verify and provide a working link or correct the citation.
Response. We verified the citation, which is correct: Bull Russ State Med Univ 2024, issue 6, pages 98-104. The DOI resolves in a browser but is refused to automated agents, which is the likely source of the difficulty.
Action. The journal archive URL was added alongside the DOI in the reference list.
- Is the codon usage of these phages similar to that of S. aureus?
Comment. Phage replication is expected to be more efficient in a host with a similar codon usage profile. Discussing this could add depth to the host-phage dynamics.
Response. This was a substantive omission and the analysis produced a result we had not anticipated. Relative synonymous codon usage was computed for each phage genome and compared with the host over the 59 codons of degenerate synonymous families. Codon usage is broadly host-like across the dataset (Pearson r 0.844-0.972), but the dataset separates completely into two groups with no intermediate value: the four genomes in which de novo prediction detects no tRNA are the four most similar to the host (r = 0.927-0.972), while the 18 tRNA-carrying genomes are the least similar (r = 0.844-0.881). The ordering is not explained by base composition - EW and SA13 depart from the host GC3 by more than seven percentage points yet resemble the host more closely in codon preference than the macro-class genomes, whose GC3 lies within four points of the host. We are explicit that the relationship is correlational and that its direction cannot be determined from sequence data.
Action. A new Methods subsection, “Codon usage and phage-encoded tRNAs”, and a new Results subsection were added; the Discussion paragraph on phage-encoded tRNAs was rewritten around the result; the per-genome data are in Supplementary Table S2.
- What constitutes the host-recognition cassette, and its domains, in the Introduction
Comment. In the Introduction, the lysis cassettes are mentioned, but what constitutes the host-recognition cassette? It would be beneficial to briefly introduce these domains in the Introduction as well.
Response. Agreed. The Introduction previously described the lysis cassette without its counterpart.
Action. A passage was added to the Introduction defining receptor-binding proteins as the tail fibres, tailspikes and baseplate-associated proteins that make first contact with the bacterial surface, noting that in Staphylococcus these recognise wall teichoic acid and its glycosylation pattern, and describing the modular N-terminal anchoring and C-terminal specificity architecture. This also answers the reviewer’s later request to define the RBP genes and to explain how the tail fibre determines host range.
- Differentiate the lysis and host-recognition modules in Table II
Comment. In Table II, it would be helpful to visually or categorically differentiate between the lysis and the host-recognition modules.
Action. The three host-recognition columns of Table II are now shaded, separating them from the lysis columns.
- Does “intact” refer to length or to 100% identity?
Comment. Lines 31-32 (Page 9): does “intact” refer to the protein’s length or to 100% sequence identity? Please specify.
Response. Neither: it refers to the completeness of the domain architecture in the mature protein.
Action. A definition sentence was added to Methods: an endolysin is described as intact when its catalytic and cell-wall-binding domain architecture is complete in the mature protein, whether or not that protein is encoded by a single uninterrupted open reading frame; the term refers to domain architecture and not to sequence identity, which is reported separately in Supplementary Table S4.
Move a portion of the descriptive lysin detail to Supplementary Material
Comment. To improve the manuscript’s readability and dynamism, I would suggest
- moving a portion of these descriptive details to a Supplementary Material.
Comment. To improve the manuscript’s readability and dynamism, I would suggest moving a portion of these descriptive details to a Supplementary Material.
Action. The two case-by-case Results paragraphs were condensed into one and the detail moved to Supplementary Text S1. All seven cases remain named in the main text with their mechanisms; only the per-case amino-acid evidence was moved. The Results section is shorter by 278 words with no loss of finding.
- A more specific comparison of vB_SauP-436A and vB_SauM-515A1
Comment. Elaborating on this comparison would strengthen the study’s narrative connection to the introduction.
Response. The two phages sit at opposite extremes of the codon usage distribution: vB_SauP-436A1 is the most host-like genome in the dataset (r = 0.972) and encodes no tRNA, whereas vB_SauM-515A1 is the least host-like (r = 0.844) and carries the three-tRNA macro-class complement. They also differ more than eight-fold in coding capacity and in endolysin architecture, the latter being intron-split in vB_SauM-515A1. We present these as genomic axes on which the two differ, not as an explanation of the biofilm phenotype, which a sequence-based analysis cannot supply.
Action. The comparison was added to the Discussion and the corresponding conclusion retained in the Abstract.
- References 24 and 26 do not seem appropriate
Comment. Line 46 (Page 9): References 24 and 26 do not seem appropriate to support this statement.
Response. Agreed on both counts. Neither reference supported the sentence as written.
Action. The citation was removed from that sentence. Reference 26 is now cited in the Introduction, where it supports the account of host-range redirection through tail-fibre modification, and reference 30 was relocated to the Discussion sentence on modular exchange of binding specificity. Neither reference is orphaned.
- Would an intron-split endolysin have the same catalytic efficiency?
Comment. Could there be structural or kinetic differences influencing the biofilm response? If the efficiency is indeed the same, then the biofilm’s behaviour must be attributed to other factors.
Response. We cannot settle this, and we say so. Splicing restores the CHAP-Amidase-SH3b architecture at the RNA level and the mature protein is reported as full length, so no structural basis for reduced activity is apparent; but neither specific activity nor turnover has been measured for the spliced product, and any residual difference would also depend on splicing efficiency in vivo, which is likewise unmeasured. The consequence for our interpretation is the same either way, as the reviewer notes: if the enzyme is equivalent, the atypical biofilm response of vB_SauM-515A1 must be attributed to factors outside the lysis module; if it is not, the difference is a kinetic property of the enzyme rather than the absence or truncation of a module, and would still be invisible to a genomic survey.
Action. The passage above was added to the Discussion, together with a sentence in the limitations stating that where the manuscript reports an intact domain architecture it reports genomic capacity and not catalytic equivalence. We name the experiment that would settle it: purified-enzyme kinetics on the spliced product.
- peer review recommendation: accept
REVIEWERS’ COMMENTS
REVIEWER #1
Here are a few minor considerations I still have:
-Species names should be italicized throughout the text.
-Did the codon usage analysis use only one S. aureus genome as a reference? I understand that the genome of the phage host strain is not always available, but its use would allow for a more accurate analysis.
-The title of Table II contains explanatory text that should be moved to footnotes.
-The citation order of the Supplementary Tables appears to be incorrect. For instance, ST1 is cited on page 22 of the proof, whereas ST3 and ST4 are first cited on pages 17 and 19, respectively.
-On page 26 of the proof, the sentence “InterPro domain support was recovered for the RBP candidates of 11 of 22 genomes” is repeated twice in a row.
-In the conclusions, the statement “Endolysin deficiency does not explain the biofilm paradox of vB_SauM-515A1” appears without prior context, making it seem disconnected.
-Finally, I still find the manuscript overly dense and detailed. The main text could be more concise, with finer details moved to the supplementary material. Furthermore, when tables are cited, the text exhaustively describes their data rather than simply highlighting the most relevant findings. While I believe the writing could still be streamlined, I do not consider this an impediment to accepting the work.
REVIEWER #2
Thank you for the revision. You’ve correctly addressed all four major comments.
Objectives/Results/Discussion/Conclusion are consistent throughout. No further round needed from me.
One minor fix before production: a paragraph in Results (RBP InterPro support, 11/22 per arm) is duplicated almost verbatim.
Manuscript is ready from my side once that’s fixed.
AUTHORS’ RESPONSE TO THE REVIEWERS
Dear Dr Vicente,
We thank both reviewers for a careful second reading. Every point raised has been addressed. One of them — the duplicated sentence in the Results, which both reviewers identified independently — was an editing error on our part that survived the first revision, and we are grateful it was caught before production.
All changes in this revision are highlighted in the manuscript file. Highlighting from the first revision was removed, so everything marked in the present version is new to this round.
Yours sincerely,
Ivan Armadi Hasugian, on behalf of both authors
Summary of changes in this revision – 2
Species names — 99 occurrences italicised across text, tables, footnotes, figure legend and supplementary captions.
Methods, codon usage — Now names the five pooled S. aureus reference genomes instead of referring to a single reference genome.
Table II caption — Reduced to a title; explanatory sentence moved to new footnote g.
Supplementary tables — Renumbered by order of first citation (VIRIDIC to S1, domain identity to S2, Phynteny to S3, codon usage to S4); caption blocks physically reordered; all nine in-text citations updated.
Results — Duplicated sentence and duplicated trailing clause removed from the annotation-source subsection.
Abstract — Orphan sentence on vB_SauM-515A1 removed from MAIN CONCLUSIONS; abstract now 190 words.
Unchanged — No figure, table value, reference or scientific claim was altered in this revision.
Reviewer 2
We address Reviewer 2 first because the single point raised is the same one Reviewer 1 also identified.
Duplicated paragraph in Results
Comment. One minor fix before production: a paragraph in Results (RBP InterPro support, 11/22 per arm) is duplicated almost verbatim.
Response. Confirmed and corrected. The defect was narrower than a whole paragraph but had the same effect: in the annotation-source subsection, the sentence beginning “InterPro domain support was recovered for the RBP candidates of 11 of 22 genomes” appeared twice consecutively, and the clause that followed it — “lower than the keyword-based call in every case, consistent with the sparse representation of receptor-binding protein folds in profile-HMM databases...” — was also repeated. Both duplications arose when two versions of the sentence were merged during the first revision.
Action. The sentence and the trailing clause now appear once. No figure, number or claim changed; only the repetition was removed. Every reported value in that passage (11 of 22 per arm, eight supported in all three arms, union 14 of 22) is unaltered.
We also note Reviewer 2’s confirmation that the Objectives, Results, Discussion and Conclusion are now consistent, and that no further round is required from their side. We have made no other change to the substance of the manuscript in this revision.
Reviewer 1
1. Species names should be italicised throughout the text
Comment. Species names should be italicized throughout the text.
Response. Corrected throughout. Italicisation was applied to every occurrence of the genus and species names in the running text, table cells, table footnotes, figure legend and supplementary captions, covering Staphylococcus aureus and its abbreviated form, Staphylococcus in “Staphylococcus phage”, and the three genera cited in the cross-genus passage of the Discussion (Shewanella oneidensis, Pseudomonas aeruginosa, Pseudomonas fluorescens).
Action. Ninety-nine occurrences were italicised. Genome accession identifiers and phage strain names were deliberately left in roman type, as these are not taxonomic names.
2. Did the codon usage analysis use only one S. aureus genome as a reference?
Comment. Did the codon usage analysis use only one S. aureus genome as a reference? I understand that the genome of the phage host strain is not always available, but its use would allow for a more accurate analysis.
Response. It did not, and we are grateful for the question because it exposed an inaccuracy in our Methods rather than in the analysis. The host reference codon table was computed by pooling the coding sequences of five curated S. aureus reference genomes — NCTC 8325 (NC_007795.1), Mu50 (NC_002758.2), MRSA252 (NC_002952.2), USA300_FPR3757 (NC_007793.1) and CP023390.1 — not from a single strain. Our Methods text referred to “the S. aureus reference genome” in the singular, which understated what was done and reasonably led to the reviewer’s question.
Response. We pooled deliberately rather than using one genome, for the reason the reviewer identifies: the propagating strain is not reported for most of these phage deposits, so no single strain can be claimed as the correct reference. Pooling across five strains that span the diversity of the species (methicillin-sensitive and methicillin-resistant, hospital- and community-associated lineages) gives a species-level codon table rather than a strain-level one, which is the level at which our comparison is made. We agree that using each phage’s own propagating host would be more accurate where that strain is known and sequenced, and this is stated as a limitation.
Action. The Methods sentence now names the five reference genomes and states explicitly that they were pooled rather than used singly.
3. The title of Table II contains explanatory text that should be moved to footnotes
Comment. The title of Table II contains explanatory text that should be moved to footnotes.
Response. Agreed. The caption carried a second sentence explaining which column rested on which kind of evidence. That sentence also duplicated material already present in footnotes e and f, so moving it removes a redundancy as well as shortening the title.
Action. The Table II caption is now a title only: “Domain architecture of lysis modules and presence of host-recognition modules among the 22 bacteriophages.” The explanatory sentence has been moved to a new footnote g.
4. The citation order of the Supplementary Tables is incorrect
Comment. The citation order of the Supplementary Tables appears to be incorrect. For instance, ST1 is cited on page 22 of the proof, whereas ST3 and ST4 are first cited on pages 17 and 19, respectively.
Response. Correct, and the problem was systematic rather than confined to those two tables. The four supplementary tables were numbered in the order in which they were written during the first revision, not in the order in which the text reaches them. First citation ran S3, S4, S1, S2.
Response. We renumbered by order of first citation rather than moving passages, because the positions of the citing passages are themselves determined by the argument: the intergenomic-similarity matrix is first needed in the dataset-curation paragraph of the Methods, and the domain-identity matrix in the paragraph defining “intact”. Reordering the text to fit the old numbering would have damaged the argument to preserve an arbitrary label.
Action. The supplementary tables were renumbered as follows, and the caption blocks were physically reordered to match: intergenomic similarity (VIRIDIC) is now Supplementary Table S1; endolysin per-domain pairwise identity is now S2; Phynteny residual annotation is now S3; codon usage similarity is now S4. All nine in-text citations were updated. First citation now runs S1, S2, S3, S4 in order. Supplementary Text S1 is a separate series and is unchanged.
5. A sentence in the Results is repeated twice in a row
Comment. On page 26 of the proof, the sentence “InterPro domain support was recovered for the RBP candidates of 11 of 22 genomes” is repeated twice in a row.
Response. Confirmed and corrected. This is the same defect Reviewer 2 identified; our response is given in full above, under Reviewer 2.
Action. The repetition and an associated repeated clause were removed. No reported value changed.
6. The Conclusions sentence on vB_SauM-515A1 appears without prior context
Comment. In the conclusions, the statement “Endolysin deficiency does not explain the biofilm paradox of vB_SauM-515A1” appears without prior context, making it seem disconnected.
Response. Agreed. The sentence sat in the MAIN CONCLUSIONS section of the structured abstract, where neither vB_SauM-515A1 nor the biofilm paradox had been introduced — the abstract reaches that sentence with no antecedent. It does not appear in the Conclusion section of the main text.
Response. We considered adding the missing context, but the abstract was already at the 200-word limit, so context could only be bought by removing a primary finding. The statement is a secondary result that is fully developed with its context in the Discussion, where the comparison between vB_SauP-436A1 and vB_SauM-515A1 is set out. We therefore removed it from the abstract rather than compress a primary finding to accommodate it.
Action. The sentence was deleted from MAIN CONCLUSIONS. The abstract is now 190 words. The finding itself is unchanged and remains in the Discussion with its full context.
7. The manuscript remains dense; details could move to supplementary material
Comment. Finally, I still find the manuscript overly dense and detailed. The main text could be more concise, with finer details moved to the supplementary material. Furthermore, when tables are cited, the text exhaustively describes their data rather than simply highlighting the most relevant findings. While I believe the writing could still be streamlined, I do not consider this an impediment to accepting the work.
Response. We accept the criticism and have not made further cuts in this revision, for a reason we would like to state rather than leave implicit. In the first revision we moved two case-by-case Results paragraphs to Supplementary Text S1, reducing that section by 278 words at the reviewer’s earlier request. The remaining density is concentrated in passages that exist because of the first round: the three-arm comparison, the domain-partitioned identity result and the intergenomic-similarity analysis were each added in direct response to specific comments, and each carries the numbers that support a claim now central to the paper.
Response. Cutting further would mean removing either those numbers or the qualifications attached to them, and our concern is that a shorter text would reintroduce exactly the kind of unsupported general statement that the first review asked us to eliminate. Given that the reviewer does not consider this an impediment to acceptance, we have preferred to leave the evidence in place. We would of course act on specific passages if the editor identifies any as expendable.
Action. No change. We note the point and would welcome specific guidance at the proof stage if further condensation is wanted.
- peer review recommendation: accept
REVIEWERS’ COMMENTS
REVIEWER #1
All my questions were addressed.
REVIEWER #2
No other comments.
- peer review recommendation: accept
