Defense-system syndromes and prophage-associated anti-phage arms races

This concept examines whether bacterial anti-phage defense systems form recurring combinations and whether broader defense repertoires are associated with prophage burden. The evidence comes from the Pan-Bacterial Anti-Phage Defense Arsenal, which analyzed seven defense-system families across the KBase Data Lakehouse pangenome. The Prophage-AMR Co-mobilization Atlas extends this question to associations between prophage markers and antimicrobial-resistance (AMR) repertoire breadth, while Prophage Gene Modules and Terminase-Defined Lineages Across Bacterial Phylogeny and Environmental Gradients tests prophage-module composition across bacterial phylogeny, environmental gradients, and NMDC metagenomes. The SNIPE Defense System in the KBase Data Lakehouse Pangenome adds the SNIPE defense system, its proposed ManYZ-dependent phage-entry mechanism, and evidence for a widespread, predominantly accessory defense family. ₁ The COG Functional Category Analysis contributes a composite functional-category signal for genes annotated as both mobile and defense-related, and discoveries records cross-project SNIPE observations, including an archaeal homologue with fitness data. _2,3

Pan-bacterial defense repertoire

Markers for at least one defense system were detected in 27,626 of 27,690 species-level pangenomes (99.8%), yielding 930,573 marker hits across 825,476 unique gene clusters. ₁ The quality-filtered analysis included 7,323 species with at least five genomes. ₁ This at-least-five-genome analysis set is the one used for the arms-race and syndrome tests, because reliable core and accessory calls require at least five sequenced genomes per species. ₁

Species-level prevalence, computed across all 27,626 pangenome species with at least one defense hit, ranged from 7.2% for cbass to 96.1% for crispr-cas; the other systems were BREX at 80.1%, restriction-modification (R-M) Type I at 76.8%, DISARM at 58.4% under the strict call, Retron at 54.7%, R-M Type II at 38.9%, and Gabija at 22.8%. ₁ Restricted to the nine phyla with at least 100 species in the at-least-five-genome analysis set of 7,323 species, the phylum-level patterns were consistent with the pan-bacterial totals: R-M Type I exceeded 80% in Pseudomonadota and Bacillota, CBASS remained below 20% in all phyla except Bacteroidota at approximately 10–15%, and Gabija reached 30–40% in Bacillota versus below 10% in Actinomycetota. ₁

The apparent CRISPR-Cas prevalence depends strongly on detection method: EggNOG description matching produced 96% species prevalence, whereas the specific Cas1 Pfam marker PF01867 produced approximately 55% prevalence on the same pangenome. ₁ The report therefore treats 96% as an upper bound, and syndrome comparisons that use the combined EggNOG-or-Pfam presence set should be interpreted as recall-oriented rather than maximally specific. ₁

SNIPE provides a complementary example of broad but annotation-sensitive defense distribution. Its DUF4041/PF13250 marker occurred in 4,572 gene clusters across 1,696 species and 33 bacterial and archaeal phyla; 86.7% of those clusters were accessory or singleton (13.3% core, 30.7% accessory, and 56.1% singleton). ₄ This supports the broader interpretation that defense markers span diverse taxa while often residing outside stable core inheritance, but it also refines prevalence estimates because only 54 of the 4,572 DUF4041 clusters carried the Mug113 description associated with the PF13455 nuclease domain. ₄ The project report presents the 4,572 clusters across 1,696 species and 33 bacterial and archaeal phyla as substantially expanding the >500 homologues reported in the original paper, with 86.7% of hits accessory or singleton. ₄

The domain identifications behind these counts are explicit. DUF4041 corresponds to Pfam PF13250 (pf13250 duf4041) and InterPro IPR025280, renamed “SNIPE associated domain”, for which the project’s domain table lists 1,612 proteins and 2,466 species; the SNIPE nuclease is PF13455 (pf13455 mug113, Mug113), which sits in GIY-YIG clan CL0418 but is not the canonical GIY-YIG family; and the E. coli SNIPE protein A0A0A1A5Z2 (558 aa) is annotated by InterPro with PF13250 at positions 232–333 and PF13455 at positions 443–520. ₄

A related annotation-specific caveat applies to prophage burden. In a separate analysis, prophage-associated annotations occurred in all 27,702 analyzed bacterial species, but the calls were based on eggNOG functional annotations rather than dedicated prophage detection; near-universal prevalence may therefore include domesticated remnants and bacterial homologs of phage genes. ₅ This refines interpretation of the defense–prophage associations below: the burden signal is robust as a marker-defined association, but true intact-prophage prevalence remains uncertain. ₅

The scale of that annotation query is documented: querying 93M eggNOG annotations with 112 WHERE clause conditions identified 4,005,537 prophage-associated gene clusters across 27,702 species. ₅

The report itself contrasts its finding that all 27,702 species carry prophage annotations with Touchon et al.’s (2016) estimate that ~46% of bacterial genomes harbor prophages, attributing the discrepancy to its broader definition: it identifies individual prophage-associated gene clusters, including domesticated remnants, not intact prophages. ₅

Defense-system syndromes

A syndrome is a recurring combination of defense systems across species. Using 1,000 phylum-stratified column-permutation null permutations and Benjamini–Hochberg false-discovery-rate (BH-FDR) correction, 27 of 28 defense-system pairs showed significant positive co-occurrence at q < 0.05. ₁ This broad positive structure supports cofitness-network-architecture by showing that defense systems can form coordinated network-level architectures rather than appearing as independent traits. ₁

The strongest pair was R-M Type II × Gabija, with 2,429 observed co-occurrences versus a null mean of 1,555, z = 46.1, and odds ratio (OR) = 24.0. ₁ The report presents this pairing as a novel pan-bacterial syndrome at this scale: R-M Type II uses sequence-specific cleavage, whereas Gabija uses a nuclease/helicase pair triggered by nucleotide-pool depletion, and both systems act on double-stranded DNA. ₁ Expressed as species fractions, the two systems co-occur in 2,429 of 4,096 R-M Type II-carrying species (59 %) versus an expected 1,555 (38 %). ₁ The report qualifies the novelty claim for R-M Type II × Gabija (OR = 24, z = 46) as being to its knowledge, cites Cheng et al. (2021, Sci Adv 7:eabe5470) for the Gabija nuclease/helicase mechanism, and proposes that the pairing may reflect complementary layered defense against different phage life-cycle stages — a possibility it offers rather than tests. ₁

Ranked by z-score, the nine leading pairs were R-M Type II × Gabija with 2,429 observed co-occurrences, null mean = 1,555, z = 46.1, and OR = 24.0; then BREX × DISARM with 4,968 observed co-occurrences, null mean = 4,629, z = 31.1, and OR = 8.2; BREX × Retron with 5,021, 4,681, z = 29.7, and OR = 8.6; DISARM × Retron with 4,280, 3,844, z = 28.4, and OR = 5.2; R-M Type I × DISARM with 4,970, 4,673, z = 27.5, and OR = 8.4; R-M Type I × BREX with 5,963, 5,732, z = 26.9, and OR = 8.7; R-M Type II × DISARM with 3,424, 2,981, z = 24.4, and OR = 4.1; R-M Type II × BREX with 3,929, 3,637, z = 23.0, and OR = 6.4; and R-M Type I × R-M Type II with 3,948, 3,665, z = 22.5, and OR = 7.6. ₁

SNIPE refines this syndrome framework by linking defense to a host entry route rather than simply adding another co-occurrence marker. Saxton et al. (2026) showed that SNIPE constitutively localizes to the inner membrane and cleaves phage DNA as it passes through the ManYZ mannose transporter pore, potentially retaining mannose transport while blocking phage infection. ₄ In E. coli K-12, direct RB-TnSeq (random barcode transposon sequencing) data from 168 experiments per ManXYZ gene show the largest defects on D-glucosamine and D-mannose, while cofitness correlations of 0.851 for manX↔manZ, 0.705 for manX↔manY, and 0.725 for manY↔manZ support their operation as a single operon. ₄ This supports a defense–metabolic trade-off hypothesis, but the transporter mechanism is established directly only for the SNIPE/ManYZ model and should not be generalized to the other syndromes. ₄

Published knockout genetics give that trade-off its quantitative shape. In E. coli, the manY (EIIC) knockout JW1807 showed zero phage growth after 24h together with loss of mannose transport and a partial glucose transport defect, and the manZ (EIID) knockout JW1808 showed a 5 orders of magnitude reduction in phage growth with the same transport cost (Burmeister et al. 2021). ₄ In phage–bacteria coevolution experiments, man mutants sweep to >95% frequency under phage pressure but decline sharply when phages evolve alternative injection mechanisms (Burmeister et al. 2021). ₄

The report’s own reading is that SNIPE resolves this resistance-versus-metabolic-cost trade-off: rather than losing ManYZ entirely, SNIPE-bearing strains retain full transporter function while cleaving phage DNA during injection, which would supply the phage-resistance benefit of man knockouts without the metabolic cost. ₄ That passage reports no direct comparative fitness measurement between SNIPE-bearing and man-knockout strains, so the resolution should be read as the report’s proposal rather than a demonstrated result. ₄

The SNIPE domain correction also refines marker-based syndrome interpretation: PF13250 is the corrected Pfam assignment for DUF4041, whereas PF13455/Mug113 is the SNIPE nuclease family and is distinct from canonical GIY-YIG PF01541; zero surveyed clusters contained both DUF4041/PF13250 and PF01541. ₄ The recommended SNIPE search therefore combines PF13250, DUF4041, or T5orf172 with PF13455 rather than treating PF01541 as the nuclease marker. ₄ The cross-project digest adds that PF13455 and PF01541 are distinct families within the same GIY-YIG clan CL0418, so the absence of DUF4041 + PF01541 clusters is biological truth rather than an annotation artifact, as verified on the E. coli SNIPE protein A0A0A1A5Z2 (PF13250 at positions 232–333, PF13455 at positions 443–520). ₃

One archaeal homologue makes the SNIPE model newly testable: the kescience-fitnessbrowser contains a methanococcus-maripaludis JJ protein (locus MMJJ_RS01635) annotated with PF13250 (DUF4041), PF13455 (Mug113/GIY-YIG nuclease) and PF10544 (DUF2525), reported as a full two-domain SNIPE protein with 129 experiments of fitness data and the first SNIPE homologue with genome-wide knockout phenotypes available for analysis. _3,4

The measured phenotype of that homologue is mild: the M. maripaludis SNIPE gene shows fitness defects under specific conditions, with minimum fitness = -1.16 on formate/acetate, but is dispensable under most conditions, and as an archaeon its phage biology context differs significantly from the Enterobacterales phage lambda system. ₄

Coverage of the nuclease domain in that resource is otherwise thin: among the 48 Fitness Browser organisms, 7 genes total carry PF13455, spread across 6 organisms (Azospirillum brasilense, Bacteroides thetaiotaomicron, Cupriavidus, Desulfovibrio vulgaris, Paraburkholderia kururiensis, and M. maripaludis), and only M. maripaludis has the full two-domain PF13250 + PF13455 architecture. ₄

The phagefoundry klebsiella genome browser reports 1 DUF4041 annotation (3 proteins at 558 aa, the same length as the E. coli SNIPE protein) and 4,619 PTS_EIIC (PF02378/ManY family) annotations, making Klebsiella the only PhageFoundry species carrying both SNIPE and the ManX-family PTS domain PF00358; the accompanying suggestion that SNIPE could affect phage-therapy efficacy in this pathogen is speculative rather than demonstrated. ₃

The underlying four-species annotation table records, for Klebsiella, 1 DUF4041 description over 3 proteins of 558 aa each, 28 mannose PTS IID (ManZ) eggNOG proteins, 1 mannose PTS IIA (ManX) eggNOG protein, 4,619 PF02378/PTS_EIIC (ManY family) Pfam annotations and 571 PF00358/PTS_EIIA_1 (ManX family) Pfam annotations. ₄ The three comparator species carry 0 DUF4041 and 0 PF00358 annotations each, but abundant PTS_EIIC annotations — 891 in Acinetobacter, 1,067 in P. aeruginosa and 260 in P. viridiflava; Klebsiella was thus the only one of the four compared PhageFoundry species in which DUF4041 was detected at all, and the report states that Klebsiella is the only PhageFoundry species with both SNIPE and the PF00358 ManX-family PTS domain by exact Pfam accession. ₄

The three Klebsiella DUF4041 proteins — JJW41_17025, KFB13_RS08150 and KFB31_RS01490 — are all 558 aa, the same length as the characterized E. coli SNIPE protein A0A0A1A5Z2, so their full-length two-domain (PF13250 + PF13455) architecture is inferred from length rather than directly confirmed here. ₄

CRISPR-Cas × CBASS was the only pair that did not reach significance, with z = 0.21 and empirical p = 0.98. ₁ These systems also had the largest prevalence gap, at 96% versus 7%. ₁ The null-model results establish co-occurrence, but they do not by themselves establish biochemical cooperation, shared genomic localization, or direct protection against the same phages. ₁

Prophage-associated defense expansion

Species-level defense-system count correlated positively with prophage-cluster burden, with Spearman ρ = 0.609, p ≈ 0, and n = 7,323. ₁ After residualizing on log10 median genome size and phylum, the partial correlation remained ρ = 0.301 with p = 1.6 × 10⁻¹⁵³. ₁ A negative-binomial generalized linear model found positive contributions from prophage-cluster count, β = 2.0 × 10⁻⁴, p < 0.001, and log10 genome size, β = 0.755, p < 0.001. ₁ Each 10-fold genome-size increase was associated with e^0.755 ≈ 2.1× more defense systems. ₁

The prophage ecology analysis supports the environmental relevance of this burden association while qualifying its interpretation. In a 1,773-species subsample, PERMANOVA (permutation-based multivariate analysis of variance) on Bray–Curtis distances of prophage module composition, with 99 permutations, found significant effects of genome-size quartile, environment, and family-level phylogeny, each with p=0.01; genome size was dominant (F=212.99), followed by environment (F=30.04) and phylogeny (F=6.17). ₅ Genome size is also the dominant univariate predictor, correlating with prophage cluster count at rho=0.717. ₅ Kruskal-Wallis tests confirmed that environmental effects remained significant within every genome-size quartile, with all p < 6.5e-78, and among 2,008 species with at least five embedded genomes, AlphaEarth environmental niche breadth correlated with prophage module count after controlling for genome size (partial Spearman rho=0.468, p=8.41e-110). ₅ Thus, the existing defense–prophage relationship is consistent with ecological structure beyond genome size and host phylogeny, although the new analysis also reports genome size as the dominant predictor and cannot exclude residual confounding. ₅

The report states that limit explicitly: genome size is the dominant predictor of prophage burden (rho=0.717), and although it is controlled via stratification and partial correlations, residual confounding cannot be excluded because larger genomes have more genes of all types, not just prophage genes. ₅

The NMDC metagenome bridge supports an abiotic-gradient component of this ecological signal: taxonomy-based inference of prophage burden across 6,365 NMDC metagenomic samples, with 87.2% median matching coverage, yielded 57 significant module–abiotic correlations at FDR < 0.05. ₅

The strongest reported Spearman correlations were packaging versus pH at 0.519, and all modules versus pH at 0.474, temperature at 0.399, depth at 0.361, and total nitrogen at 0.333. ₅

The report offers an interpretation of the pH signal that it does not test: the positive correlation with pH (rho=0.519 for packaging, 0.474 overall) is read as alkaline environments possibly favouring lysogeny, consistent with stress-induced prophage induction being pH-sensitive. ₅

Cross-validation between pangenome enrichment and the NMDC correlations was concordant: B_head_morphogenesis, C_tail, and G_anti_defense were significant in both analyses. ₅

Four additional modules — A_packaging, D_lysis, E_integration and F_lysogenic_regulation — are NMDC-significant but not pangenome-enriched beyond phylogenetic expectation, which the report interprets as their near-universal presence masking environment-specific variation in the pangenome analysis. ₅

The AMR co-mobilization analysis supports the direction of this prophage-burden pattern: across 4,770 species, prophage marker density was positively associated with AMR repertoire breadth (Spearman rho=0.572, p<10^-300, n=4,770), with a partial correlation of rho=0.464 (p=1.0×10^-253) after controlling for genome count. ₆ A log-log regression gave a slope of 0.823 (SE=0.018, p<10^-300) with R²=0.30, which the report reads as a 10-fold increase in prophage density predicting a ~6.6-fold increase in AMR breadth, and genome-normalized prophage and AMR counts correlated at rho=0.608 (p<10^-300). ₆ The association was significant across Pseudomonadota (rho=0.54), Bacillota_A (rho=0.55), Bacillota (rho=0.40), Bacteroidota (rho=0.59), and Actinomycetota (rho=0.29). ₆ This is a related species-level association, not evidence that defense systems themselves cause AMR accumulation or that prophages mediate the transfer. The report itself reads the association as consistent with two non-exclusive mechanisms — prophages directly mobilizing resistance genes via specialized or generalized transduction, and species with high recombination potential (indicated by many prophages) also being better at acquiring genes from other mobile elements — and states that it does not prove phage-mediated AMR transfer, because species with open pangenomes may independently accumulate both prophages and AMR genes. ₆

External literature is cited on both sides of that interpretive gap. The report aligns its species-level finding with Rendueles et al. (2018), who showed that capsule-encoding bacteria have more prophages and more antibiotic resistance genes across >100 pangenomes, and positions its own 4,770-species analysis as a much larger-scale extension that controls for genome count rather than capsule presence. ₆ It also cites Bearson & Brunelle (2015), who demonstrated that fluoroquinolone exposure induces prophage in multidrug-resistant Salmonella, enabling phage-mediated transduction of resistance plasmids, and Fisarova et al. (2021), who showed that S. epidermidis phages transduce antimicrobial resistance plasmids at high frequency; both are external mechanistic results, while this project’s own prophage–AMR evidence remains correlative. ₆

SNIPE supports the interpretation that mobile defense systems can contribute to this broader burden pattern: 80.4% of DUF4041-containing clusters (3,675/4,572) were represented by T5orf172 or DUF4041 descriptions, and the accessory-plus-singleton fraction was 86.7%. ₄ However, the environmental association was modest rather than habitat-diagnostic: among AlphaEarth-embedded species, SNIPE-bearing species numbered n=1,069 and non-SNIPE species n=13,977; 22 of 64 dimensions differed significantly at Bonferroni-corrected p < 0.05, with the largest effect in A19 at Cohen’s d = 0.26 and p = 5.0e-14. ₄

On that basis — 22/64 AlphaEarth dimensions at Bonferroni-corrected p < 0.05 — the report describes SNIPE-bearing species as occupying statistically distinct environmental niches and rejects the null hypothesis of random SNIPE distribution for both taxonomic distribution and environmental niche association, while the niche difference itself remains a modest effect. ₄

All nine major phyla showed positive, significant partial correlations between defense count and prophage-cluster burden: Campylobacterota, n = 102, ρ = 0.530, p = 9.8×10⁻⁹; Bacillota, n = 735, ρ = 0.481, p = 8.4×10⁻⁴⁴; Bacillota_A, n = 1,134, ρ = 0.419, p = 2.0×10⁻⁴⁹; Verrucomicrobiota, n = 129, ρ = 0.369, p = 1.6×10⁻⁵; Patescibacteria, n = 191, ρ = 0.343, p = 1.2×10⁻⁶; Pseudomonadota, n = 2,172, ρ = 0.342, p = 9.6×10⁻⁶¹; Cyanobacteriota, n = 124, ρ = 0.309, p = 4.9×10⁻⁴; Bacteroidota, n = 925, ρ = 0.261, p = 6.3×10⁻¹⁶; and Actinomycetota, n = 803, ρ = 0.185, p = 1.3×10⁻⁷. ₁ The report also notes that effect strength itself varies phylogenetically, with the Actinomycetota effect approximately 3× weaker than the Campylobacterota effect, and flags the covariates that predict this effect-size heterogeneity as unexplored. ₁

Together, these associations are consistent with a prophage-associated anti-phage arms-race hypothesis, but they do not formally establish a phylogenetically independent arms race. ₁ The analysis controlled for phylum categorically rather than using a phylogenetically corrected regression, and the negative-binomial model fixed dispersion at alpha = 1.0 rather than estimating it, a conservative default under which both focal coefficients remained significant at p<0.001 but standard errors may be slightly under- or over-stated. ₁ The prophage classifier also included broad module candidates such as integrase, holin, endolysin, CI-like repressor, and tail proteins that can match non-phage bacterial genes; consequently, n_prophage_modules saturated at 7 for 35% of species, while the unbounded n_prophage_clusters count was used as the primary arms-race predictor. ₁

Prior work frames what is new here. The report attributes to Vassallo et al. (2022, Nat Microbiol 7:1568) an average of approximately 5 defense systems per bacterial genome and substantial cross-phylum variation in per-genome defense repertoires — findings that predict, but did not directly test, the arms-race hypothesis at pangenome scale. ₁ It further states that Tesson 2022 and Vassallo 2022 reported system-level distributions and pairwise co-occurrence heatmaps but did not explicitly control for phylum plus genome size in an arms-race regression, nor test syndrome significance against a phylum-preserving null with corrected q-values. ₁

The new analysis refines the arms-race interpretation by separating near-universal core-like modules from environmentally variable functions. Packaging and lysogenic-regulation modules occurred in 100.0% of species and lysis in 99.9%, whereas anti-defense occurred in 64.3%, head morphogenesis in 56.1%, and tail in 55.6%; only 34.9% carried all 7 modules; across all 27,702 species surveyed, 4,005,537 prophage gene clusters were identified via eggNOG annotations. ₅ The report attributes the near-universality of the packaging, lysis and lysogenic-regulation modules partly to domesticated prophage remnants (Bobay et al. 2014) — defective prophages under purifying selection that retain packaging, lysis and regulatory genes — and argues that the structurally variable head-morphogenesis, tail and anti-defense modules are more informative for prophage ecology because their presence indicates more complete prophage elements. ₅ In constrained permutation null models — 500 permutations preserving host family × genome size quartile strata across 18,031 species, which identified 8 significant module × environment enrichments at FDR < 0.05 — human-associated environments were enriched for tail (log2(OR)=2.21, Z=10.86), head morphogenesis (log2(OR)=1.98, Z=10.00), and anti-defense (log2(OR)=1.70, Z=8.76), while anti-defense was depleted in freshwater (log2(OR)=-0.74, Z=-4.74) and animal-associated environments (log2(OR)=-0.24, Z=-8.77). ₅ Human-clinical environments were also enriched for anti-defense (log2(OR)=0.14, Z=5.14), tail (log2(OR)=0.77, Z=3.90), and head morphogenesis (log2(OR)=0.70, Z=3.59). ₅ These associations are consistent with stronger phage-mediated selection in some human-associated settings, but they are not a direct mechanistic demonstration of coevolution. ₅

Mapped across 142 GTDB phyla, prophage module prevalence is remarkably uniform for modules A (packaging), D (lysis) and F (lysogenic regulation), all ≥99% in the top 15 phyla, while the most phylogenetically variable modules span wide ranges: B (head morphogenesis) from 10.2% (Thermoplasmatota) to 81% (human_associated species), C (tail) from 6.5% (Thermoplasmatota) to 82% (food-associated species), and G (anti-defense) from 3.1% (Patescibacteria) to 97.8% (Campylobacterota). ₅ The report’s summary states the same split directly: not all prophage functions respond equally to environment, with the ubiquitous modules A, D and F showing minimal environmental variation while the structurally variable modules B, C and G carry the environmental signal. ₅

Environment classification of 293,050 genomes via NCBI isolation_source metadata yielded 10 categories with at least 30 species, and chi-squared tests confirm significant environment × module associations for B_head (chi2=1,051, p=1.72e-220), C_tail (chi2=1,749, p≈0) and G_anti_defense (chi2=818, p=2.58e-170). ₅

Within-pangenome conservation follows the same division: modules D (lysis, 62.6% core) and A (packaging, 57.0% core) are dominated by core genes, while modules E (integration, 59.3% singleton) and B (head morphogenesis, 56.0% singleton) are dominated by singletons, which the report reads as phage structural genes diversifying rapidly while regulatory and lysis genes remain more conserved (Hendrix et al. 2000). ₅

Two external results are cited as context rather than as measurements from this corpus: Low et al. (2024) demonstrated enrichment of prophage-encoded cargo genes in human-impacted environments across 38,605 genomes, and Tesson et al. (2024) found that 89% of anti-defense systems localize in mobile genetic elements, predominantly prophages. ₅

mmseqs2 clustering of 38,085 terminase large subunit (terl) protein sequences from 11,789 species at 70% amino-acid identity (AAI) yielded 10,991 lineages with a highly skewed size distribution (median=1, max=1,094 across 869 species); the largest lineage, L_SFJM01000051.1_13, spans 1,094 members across 869 species as a cosmopolitan prophage type, 6,921 lineages (63%) are singletons representing species-specific prophage types, and threshold sensitivity analysis gave 4,001 lineages at 50% AAI and 16,283 at 80% AAI. ₅ The new evidence also refines the phylogenetic interpretation: among 824 TerL lineages with at least 5 species, 325 were classified as specialists (Shannon entropy < 1.0 or dominant environment > 80%) and 499 as generalists distributed across 3+ environments, yet no individual lineage showed significant environment-specific enrichment in 0/500 tests after FDR correction. ₅ Module-level environmental ecology and lineage-level ecology therefore appear decoupled in this analysis, consistent with modular exchange without independent whole-lineage adaptation. ₅

The report positions that null result against prior work: the lack of significant lineage-level enrichment contrasts with Paez-Espino et al.’s (2016) finding of strong habitat-type specificity for most viral groups, but its analysis used constrained permutations preserving host family composition, under which the lineage distribution is well explained by host phylogeny alone — consistent with Mavrich & Hatfull (2017), who showed that temperate phage evolution is primarily shaped by host range rather than environmental selection. ₅

Accessory acquisition and interpretation limits

Six of seven systems were significantly enriched in the auxiliary and singleton pangenome relative to a background of 46.8 % core and 37.9 % singleton gene clusters across 132,531,501 gene clusters, with χ² p ≈ 0 for all systems. ₁ The report counts six of seven system families as accessory-enriched while its per-system table carries eight rows, because R-M Types I and II are listed separately there. ₁ Retron had 8.2% core and 65.8% singleton fractions; Gabija had 9.1% and 68.8%; R-M Type II had 9.4% and 68.2%; R-M Type I had 9.8% and 65.3%; CBASS had 11.8% and 61.1%; BREX had 20.7% and 55.8%; CRISPR-Cas had 27.4% and 47.2%; and DISARM had 40.6% and 42.3%. ₁

For the accessory-enriched systems, singleton fractions were 1.5–1.8× above the 37.9% background baseline, which is consistent with horizontal transfer of defense loci through mobile genetic elements. ₁ This supports gene-function-acquisition-depth by linking defense-system presence to accessory and singleton gene-cluster distribution, while the inference remains limited by marker specificity. ₁

A functional-category signal supports the mobile-defense-island framing from an annotation angle rather than a marker angle: in the COG (clusters of orthologous groups) analysis, the composite LV (mobile+defense) category showed +0.34% enrichment with 76% consistency. _2,3 Those reports propose functional modules such as “mobile defense islands” as a possible interpretation of the LV signal rather than demonstrating such modules directly, and argue that these composite-category genes should not be filtered out as noise because they represent genuine multi-functional genes. _2,3

The AMR inventory supports the accessory-carriage picture at marker scale: across the full pangenome inventory, 83,008 AMR gene clusters and 3.47 million prophage marker clusters were identified; prophage markers were overwhelmingly accessory (83.8%) and frequently singleton (53.8%), while AMR genes were also predominantly accessory (69.7%) but less singleton (36.1%), and 14,669 of 27,702 species (52.9%) carried both AMR and prophage markers. ₆

SNIPE supports this mobile-defense interpretation: among its 4,572 DUF4041-containing clusters, 13.3% were core (604 clusters), 30.7% accessory (~1,402 clusters), and 56.1% singleton (2,566 clusters), and the marker spanned 1,696 species and 33 phyla. ₄ The report reads the combined 86.7% accessory-plus-singleton fraction as strong support for the defense island mobility hypothesis, with only 13.3% of SNIPE clusters core within their species pangenome indicating that SNIPE is typically gained or lost rather than vertically inherited; horizontal transfer is inferred rather than directly measured here, DUF4041 can occur in non-SNIPE proteins, and eggNOG annotations may miss divergent homologues. ₄

DUF4041-bearing species are spread unevenly across phyla, with 556 species in Pseudomonadota, 334 in Actinomycetota, 275 in Bacillota_A, 206 in Bacillota, 114 in Bacteroidota, 45 in Nitrospirota, 28 in Cyanobacteriota and 22 in Planctomycetota; the report reads this patchy distribution across many phyla and families as consistent with, but not proof of, horizontal gene transfer of defense islands. ₄ The cross-project digest records the same survey of DUF4041 (PF13250) in the 293K-genome KBase Data Lakehouse pangenome — 4,572 gene clusters across 1,696 species and 33 phyla, 86.7% of them accessory or singleton — and notes that this far exceeds the >500 homologues reported in Saxton et al. 2026 (Nature), with the patchy distribution across many phyla and families supporting horizontal gene transfer. ₃

The prophage ecology analysis refines the mobile-module model: packaging, lysis, and lysogenic regulation were more consistently co-localized, whereas integration and anti-defense were more dispersed across contigs, with mean co-localization values of 0.986, 0.901, and 1.000 for the three core-like modules versus 0.420 and 0.281 for integration and anti-defense. ₅ This pattern is consistent with a physically linked prophage backbone plus accessory defense islands, but it was tested in only 15 species and does not establish transfer between elements. ₅

That test used 15 phylogenetically stratified species, one per phylum, with at most 300 genomes each and 200 null permutations; 44 of 95 module–species tests (46.3%) showed significantly higher within-module gene co-occurrence than expected by chance at one-tailed z > 1.96, and mean contig co-localization across all tests was 0.769. ₅

Per module, the significant-test count, mean Z-score and mean contig co-localization were: packaging 11/15, 3.40, 0.986; head morphogenesis 1/6, -0.44, 0.474; tail 5/12, 1.43, 0.752; lysis 13/15, 4.45, 0.901; integration 1/15, -1.93, 0.420; lysogenic regulation 10/15, 2.16, 1.000; and anti-defense 3/11, 1.50, 0.281. ₅

Integration is an explicit negative result — 1 of 15 tests significant at a mean Z of -1.93 — and the report relates it to Bobay et al. (2014), who showed that domesticated prophage remnants, which are more common than intact prophages, often retain only integrases while losing the rest of the prophage backbone. ₅

The low co-localization of anti-defense genes (0.281) is aligned by the report with Pinilla-Redondo et al. (2020) and Rahimian et al. (2026), who found that anti-defense and defense systems co-localize into functional “defense islands” that are genomically distinct from the core prophage. ₅

The scope of that co-occurrence analysis is narrow: 15 species, one per phylum, give broad phylogenetic coverage but limited within-phylum replication; head morphogenesis was testable in only 6 species because of sparsity; and the 200-permutation null model provides adequate but not exhaustive precision for z-score estimation. ₅

The AMR analysis qualifies a simple immediate-neighborhood interpretation of mobile-element association. In the top-100 AMR-burdened species, 55.7% of AMR gene instances reside on contigs that also carry strict prophage markers (terminase, phage structural proteins, holin/lysin). ₆ Among 36,041 AMR gene instances, 12,026 (33.4%) were within 50 genes of a strict prophage marker, but only 3,731 (10.4%) were within 10 genes and 1,991 (5.5%) within 5 genes; the median distance among co-localized genes was 34 genes. ₆ AMR genes within 10 genes were 67.6% accessory versus 65.5% for distal genes, with OR = 1.10, one-sided p=0.005, and bootstrap 95% confidence interval [1.024, 1.185], a modest 2.1 percentage point difference. ₆ The threshold-dependent odds ratios—0.78 at 3 genes, 0.92 at 5, 1.10 at 10, 1.19 at 15, and 1.28 at 50—together with a median species-level OR of 0.85 indicate heterogeneous proximity effects rather than uniform co-mobilization. ₆ Prophage identification in that analysis used keyword and Pfam matching on bakta annotations rather than dedicated prophage prediction tools such as PHASTER or geNomad, so it may include false positives such as phage-defense systems and may miss divergent prophages. ₆

The DISARM accessory-enrichment result is unreliable because the DrmB SNF2 helicase Pfam anchor PF00176 is a widespread housekeeping-helicase domain. ₁ Retron detection likewise used the broad RVT_1 Pfam PF00078, and the stringent call reduced candidates from 15,109 to 15,098 by dropping 11 species; these candidates should therefore be interpreted as reverse-transcriptase candidates in defense-syndrome context rather than characterized retron systems. ₁

The analysis set contained 7,323 of 27,690 pangenome species, or 26%, because of the at-least-five-genome filter. ₁ This sampling restriction favors well-sampled, culturable, high-priority organisms, so extrapolation to the full 293K-genome tree, including many environmental MAGs, remains untested. ₁ The AMR analysis likewise sampled 20 genomes per species for its 100-species neighborhood analysis rather than all 293K genomes, so its proximity estimates may change under exhaustive sampling. ₆

The prophage ecology analysis adds that only 28% of genomes had AlphaEarth environmental embeddings, producing a biased subsample toward clinically and environmentally well-sampled lineages; NCBI isolation_source metadata were also sparse and inconsistently labeled. ₅ Its NMDC validation was indirect, relying on genus-level conservation of prophage content, although taxonomy-based burden inference across 6,365 samples achieved 87.2% median matching coverage. ₅ SNIPE’s AlphaEarth comparison was similarly limited to 28.4% of genomes and was biased toward environmental isolates with latitude/longitude metadata; the dimensions, including A19, have no assigned biological interpretation. ₄

Tensions

The observed breadth of the defense arsenal is partly in tension with marker-level specificity: EggNOG description matching gives 96% CRISPR-Cas prevalence, while the Cas1 PF01867 marker gives approximately 55%, and the DISARM anchor can identify non-DISARM SNF2 helicases. ₁ SNIPE adds a related tension: 4,572 DUF4041 clusters were detected, but only 54 carried the Mug113 co-annotation, and DUF4041 may occur outside complete SNIPE architectures. ₄ The co-occurrence and arms-race patterns are therefore strongest as system-presence associations, while their mechanistic interpretation requires context-aware multi-marker validation. ₁

The modest and threshold-dependent AMR–prophage proximity effect is also in tension with a strong direct-mobilization reading of the species-level prophage-density/AMR-breadth association: the latter has rho=0.572 across 4,770 species, whereas only 10.4% of AMR instances were within 10 genes of a strict prophage marker. ₆ These results need not conflict, because species-level co-acquisition and transfer need not produce immediate gene adjacency; dedicated prophage calls and sequence-resolved analyses are required to distinguish those explanations. ₆

The new prophage ecology results create a further qualification rather than a resolved contradiction: environmental effects on prophage-module composition remained after genome-size stratification and exceeded the reported family-level phylogenetic effect, but genome size was still dominant (F=212.99), and the prophage calls may include bacterial homologs and domesticated remnants. ₅ This is compatible with an ecological arms-race signal, but does not by itself establish that intact prophages caused the environmental or defense associations. ₅

Two internal overstatements in that report are recorded rather than reconciled. Its summary calls anti-defense the most environmentally responsive module, reporting human-associated enrichment and freshwater and animal-associated depletion, but its own enrichment table gives tail log2(OR)=2.21 with Z=10.86, head morphogenesis 1.98 with 10.00, and anti-defense 1.70 with 8.76, so anti-defense is not strongest by either listed measure. ₅ Its summary likewise presents the co-occurrence analysis across 15 phyla as empirical validation of prophage modularity at pangenome scale, concluding that packaging, lysis and lysogenic-regulation genes behave as physically linked functional units while integration and anti-defense genes are “genomically independent”, and supporting a two-tier model of prophage organization (core backbone plus scattered accessory functions); the stronger-linkage half of that claim matches its co-occurrence table, but the independence half overstates it, since the table reports significant co-occurrence in 1/15 and 3/11 tests with mean contig co-localization of 0.420 and 0.281 — weaker co-occurrence and more scattered placement, not demonstrated independence. ₅

SNIPE’s low Lambdavirus infection rate in PhageFoundry—1 of 188 E. coli strains, or 0.5%—is consistent with widespread ManYZ variation or loss, but this does not establish that SNIPE caused the pattern; no Klebsiella SNIPE or ManYZ fitness data were available. ₄ The report also detected DUF4041 in Klebsiella but not in Acinetobacter, P. aeruginosa, or P. viridiflava in the four-species comparison, while co-occurrence with mannose-transporter annotations was only consistent with, rather than proof of, defense at the phage-entry route. ₄

The phage in question is Lambdavirus 411_P1 (genus Lambdavirus), which infects only 1 of 188 strains (NILS06) — a 0.5% infection rate, the lowest of any phage in that dataset, which the report reads as consistent with widespread ManYZ variation or loss conferring near-universal lambda resistance rather than as a demonstrated cause. ₄ The report further attributes to Gaborieau et al. the comparative conclusion that adsorption factors (surface receptors), not antiphage defense systems, are the primary predictors of phage–host interaction, and then reads this as supporting SNIPE’s niche at the adsorption/injection step; that application is an extrapolation, since SNIPE itself was not tested in the interaction matrix. ₄

One numeric discrepancy inside the defense-arsenal report itself is recorded here rather than reconciled: its narrative section reports the Actinomycetota arms-race association as partial ρ = 0.18 with p < 10⁻⁷, noting that CRISPR-Cas is very rare in Mycobacterium and Streptomyces gene-content data, whereas its per-phylum results table reports n_species = 803, partial ρ = 0.185 and p = 1.3×10⁻⁷ for the same phylum. ₁ The narrative figure refers to the phylum-level defense-count/prophage-burden association rather than to a correlation measured specifically in Mycobacterium or Streptomyces; the table values are the more precise of the two, the p-value statements conflict, and both are kept under separate attribution instead of being averaged. ₁

Open Directions

• Apply PADLOC MacSyFinder-style multi-PFam and gene-order rules to Retron, DISARM, and Gabija, and test whether the R-M Type II × Gabija association remains at 2,429 observed co-occurrences versus a null mean of 1,555 after marker refinement. ₁
• Extend the current focused 7-system panel to 20+ defense systems using DefenseFinder/PADLOC HMMs — CRISPR-Cas subtypes, Zorya, Thoeris, Wadjet, Druantia, pAgo, PARIS, ThsA-ThsB, and defense-associated antitoxin cassettes — which would enable a full defense-repertoire principal-component analysis and genus-level clustering of “defensive phenotype”, and ask whether the 27-of-28 positive-pair pattern persists. ₁
• Fit a GTDB-tree phylogenetic mixed-effects arms-race model using defense count, n_prophage_clusters, genome size, and phylum, and test whether the residual association remains after explicit phylogenetic correction. ₁
• Apply geNomad or PHASTER to the KBase Data Lakehouse genomes and use scaffold sequences for base-pair-resolution distances, then test whether the species-level AMR breadth association persists after separating phage, plasmid, and ICE mobilization. ₆
• Re-call prophages with dedicated tools and test whether the environmental enrichment of tail, head morphogenesis, and anti-defense remains after excluding bacterial homologs and domesticated remnants. ₅
• Use sequence-resolved genomes and the 6,365-sample NMDC bridge to test whether pH-associated prophage burden reflects intact prophage abundance, genus-level inference, or pH-linked community composition. ₅
• Identify the species carrying only R-M Type II versus those carrying both R-M Type II and Gabija (OR = 24), then check in phage-challenge experiments in a tractable host such as E. coli or B. subtilis whether the Gabija addition provides a quantitative fitness advantage — an advantage this corpus proposes but does not establish. ₁
• Cross-reference defense syndromes and prophage burden with NCBI isolation-source metadata or AlphaEarth embeddings to test whether habitat-specific defense architectures explain variation not captured by phylum. ₁
• Revisit prophage-proximal AMR fitness costs as fitness-browser coverage expands, including all available genomes from Klebsiella pneumoniae, Acinetobacter baumannii, Pseudomonas aeruginosa, and Escherichia coli. ₆
• Use PF13250/DUF4041 plus PF13455 gene-order and sequence-profile rules to distinguish complete SNIPE architectures from divergent or non-SNIPE DUF4041 proteins, then test whether SNIPE prevalence and accessory status remain at 4,572 clusters and 86.7% accessory-or-singleton after refinement. ₄
• Generate ManYZ and SNIPE mutant fitness and phage-challenge measurements in Klebsiella and other phage-therapy targets, testing whether transporter preservation produces the predicted resistance-versus-metabolic-cost trade-off; the report notes this would require new mutant libraries in natural SNIPE-carrying strains such as K. pneumoniae NCTC9140 or HS11286. ₄
• Use NMDC longitudinal samples, if available, to test whether prophage burden changes over time within the same site, since the reported environmental correlations are cross-sectional and within-site temporal change remains untested. ₅