Pan-Bacterial Anti-Phage Defense Arsenal

Raw report: phage_defense_arsenal Feeds into: cofitness-network-architecture, cross-tenant-data-bridging, functional-marker-validation, gene-function-acquisition-depth, homology-search-negative-evidence, horizontal-gene-transfer-driven-innovation, pangenome-integration, phage-defense-syndromes-and-arms-race, phylogenetic-confounding-of-pangenome-associations, sampling-depth-and-downsampling-effects, scale-dependent-mobile-element-associations, two-speed-bacterial-genome

Overview

This report analyzes seven anti-phage defense-system families across the KBase Data Lakehouse pangenome. Markers were detected in 27,626 of 27,690 species-level pangenomes (99.8%), producing 930,573 marker hits across 825,476 unique gene clusters. The ≥5-genome quality-filtered analysis set contains 7,323 species and is used for the arms-race and syndrome tests. The report uses it because reliable core/accessory calls require ≥5 sequenced genomes per species. ₁

The report separates three call tiers. 27,626 / 27,690 species (99.8 %) have at least one defense hit. 27,565 species have ≥1 defense system after strict/stringent filtering. 7,323 species pass the ≥5-genome pangenome-quality threshold that defines the arms-race and syndrome analysis set. The species-level presence/absence matrix data/species_defense_matrix.tsv.gz covers 27,626 species × 8 system flags (7 systems plus Retron_candidate versus Retron_stringent), with covariates phylum, median_genome_size, no_genomes and n_clusters. ₁

The report tested 19 markers across 7 systems against both interproscan_domains and bakta_pfam_domains. Broad-Pfam markers needed anchor-plus-context filtering at species level: Retron RVT_1, DISARM PLD/SNF2 and Gabija UvrD. Specific markers can be used as single-marker calls: CRISPR Cas1, CBASS SAVED, BREX PglZ and Gabija OLD_TOPRIM_C. ₁

The report cites two resources for modern consolidated defense-system inventories and detection rules: DefenseFinder (Tesson et al. 2022, Nat Commun 13:2561) and PADLOC (Payne et al. 2021, NAR 49:10868). It also cites Vassallo et al. (2022, Nat Microbiol 7:1568), which reported that bacterial genomes carry ~5 defense systems on average and that per-genome defense repertoires vary substantially across phyla. Those findings predict the arms-race hypothesis at pangenome scale but did not directly test it. The authors say that prior work (Tesson 2022, Vassallo 2022) reported system-level distributions and pairwise co-occurrence heatmaps. They say that prior work did not explicitly control for phylum plus genome size in an arms-race regression. They also say it did not test syndrome significance against a phylum-preserving null with corrected q-values. The authors contrast their own analysis with that work on both points. ₁

Key Findings

Defense-system prevalence

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, with 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%. The report restricted phylum comparisons to the 9 phyla with ≥100 species in the ≥5-genome analysis set (7,323 species), and 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. Prevalence by phylum is plotted in figures/system_prevalence_by_phylum.png. ₁

Defense count and prophage burden

Species-level defense-system count correlated with prophage-cluster burden at Spearman ρ = 0.609 (p ≈ 0, n = 7,323). After residualizing on log10 median genome size and phylum, partial ρ remained 0.301 (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 report calls the partial correlation a moderate effect that survives the two strongest known confounders. It remains an association with prophage-cluster burden, not a direct measurement of phage pressure. figures/arms_race_scatter.png plots defense-system count against prophage-cluster burden, coloured by phylum, with an overall ρ inset. ₁

Arms-race scatter — defense count vs prophage cluster burden by phylum. From phage_defense_arsenal.

All nine major phyla showed positive, significant partial correlations: 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 notes that effect strength varies phylogenetically. The Actinomycetota effect is ~3× weaker than the Campylobacterota effect, and the report flags the covariates that predict this effect-size heterogeneity as worth exploring. figures/partial_correlation_barplot.png plots the per-phylum partial Spearman ρ with significance stars. The report’s discussion frames the result as arms-race universality across 9 phyla and says the effect is not lineage-restricted. As support, it states that Actinomycetota, where CRISPR-Cas is very rare in Mycobacterium and Streptomyces gene-content data, still shows partial ρ = 0.18 (p < 10⁻⁷). That narrative value describes the phylum-level Actinomycetota association, not a correlation measured within Mycobacterium or Streptomyces. Its p-value statement conflicts with the results table, which gives Actinomycetota n_species = 803, partial ρ = 0.185 and p = 1.3×10⁻⁷. The report does not reconcile the two, so both values are kept here as reported. ₁

Defense-system syndromes

Using 1,000 phylum-stratified column-permutation null permutations and Benjamini–Hochberg false-discovery-rate correction (BH-FDR), 27 of 28 defense-system pairs showed significant positive co-occurrence at q < 0.05. 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. Other strong pairs included BREX × DISARM (n = 4,968, null mean = 4,629, z = 31.1, OR = 8.2), BREX × Retron (5,021, 4,681, z = 29.7, OR = 8.6), DISARM × Retron (4,280, 3,844, z = 28.4, OR = 5.2), R-M Type I × DISARM (4,970, 4,673, z = 27.5, OR = 8.4), R-M Type I × BREX (5,963, 5,732, z = 26.9, OR = 8.7), R-M Type II × DISARM (3,424, 2,981, z = 24.4, OR = 4.1), R-M Type II × BREX (3,929, 3,637, z = 23.0, OR = 6.4), and R-M Type I × R-M Type II (3,948, 3,665, z = 22.5, OR = 7.6). Two figures show the co-occurrence structure. The first is a matrix of signed log₁₀(OR) for all pairs with BH-FDR significance stars (figures/syndrome_heatmap.png). The second is a circular network of positive syndromes at q < 0.05, with edge width proportional to z-score (figures/syndrome_network.png). ₁

CRISPR-Cas × CBASS was the only pair that did not reach significance (z = 0.21, p_emp = 0.98). CRISPR-Cas and CBASS also had the largest prevalence gap, at 96% versus 7%, and the report describes CBASS as essentially independent of CRISPR-Cas presence. R-M Type II and Gabija co-occur in 2,429 of 4,096 R-M Type II-carrying species (59 %), versus an expected 1,555 (38 %). The R-M Type II × Gabija pairing is the report’s strongest syndrome (OR = 24, z = 46). The authors describe it as a novel pan-bacterial finding at this scale, qualified as being to their knowledge; its two systems both act on double-stranded DNA, with Type II R-M using sequence-specific cleavage and Gabija using a nuclease/helicase pair triggered by nucleotide-pool depletion. This mechanistic rationale is cited from Cheng et al. 2021 (Sci Adv 7:eabe5470); the report does not test it. The report also proposes, as a possibility rather than an established result, that the pairing may reflect complementary layered defense against different phage life-cycle stages. ₁

Defense syndrome heatmap — signed log₁₀(odds ratio). From phage_defense_arsenal.

Accessory-pangenome enrichment

Six of seven systems were significantly enriched in the auxiliary and singleton pangenome relative to a pangenome-wide background of 46.8% core and 37.9% singleton across 132,531,501 gene clusters, with χ² p ≈ 0 for all systems. Core and singleton fractions were: Retron, 8.2% and 65.8%; Gabija, 9.1% and 68.8%; R-M Type II, 9.4% and 68.2%; R-M Type I, 9.8% and 65.3%; CBASS, 11.8% and 61.1%; BREX, 20.7% and 55.8%; CRISPR-Cas, 27.4% and 47.2%; and DISARM, 40.6% and 42.3%. The report labels five systems as strong accessory: Retron, Gabija, R-M Type II, R-M Type I and CBASS. It labels BREX as moderate accessory, CRISPR-Cas as mild accessory and DISARM as near baseline. The per-system table lists R-M Types I and II as separate rows within the seven system families. figures/core_vs_accessory_by_system.png shows stacked bars of core, non-singleton auxiliary and singleton fractions for each system against the pangenome background. ₁

For the accessory-enriched systems, singleton fractions were 1.5–1.8× above the 37.9% background baseline, consistent with active horizontal transfer of defense loci through mobile genetic elements, though not a direct demonstration of it. Such transfer is a core prediction of the “defense islands” hypothesis of Rocha & Bikard 2022. DISARM was near baseline and is treated as an unreliable accessory-enrichment result because the DrmB SNF2 helicase Pfam anchor (PF00176) is a widespread housekeeping-helicase domain. The report calls this almost certainly a detection artefact that inflates the DISARM cluster set with non-mobile SNF2 helicases. ₁

Detection and data-engineering findings

EggNOG description-based CRISPR detection produced 96% species prevalence, whereas detection using the specific Cas1 Pfam marker PF01867 produced approximately 55% prevalence on the same pangenome, a difference of approximately 40 percentage points. The report therefore treats the 96% value as an upper bound and warns that CRISPR comparisons must specify their detection method. ₁

The primary Pfam-domain source was the 833-million-row interproscan_domains table rather than the 18.8-million-row bakta_pfam_domains table. Filtering InterProScan for analysis = 'Pfam' and a curated signature list of approximately 20 accessions, followed by joining to gene_cluster, returned ~500K rows in under 90 seconds; bakta_pfam_domains contained 0 Cas1 PF01867 hits compared with ~25K in InterProScan. InterProScan uses version-free accessions such as PF01867, whereas Bakta uses versioned accessions such as PF01867.29, so cross-database joins must strip Bakta versions. The source collection was kbase_ke_pangenome (kbase-ke-pangenome). The report drew on its interproscan_domains, eggnog_mapper_annotations, bakta_pfam_domains, gene_cluster, pangenome, gtdb_metadata, gtdb_taxonomy_r214v1 and genome tables. These supplied defense-marker Pfam hits, core/auxiliary/singleton flags, phylum and genome-size covariates, and species-level pangenome statistics. ₁

Spark Connect .write.parquet() writes to cluster storage (S3) rather than the client notebook filesystem. A call like df.write.mode("overwrite").parquet("../data/foo.parquet") returns success but leaves only an empty _SUCCESS marker in the local ../data/foo.parquet/ directory. The report recommends converting results with .toPandas() and then using pandas.to_csv(..., compression="gzip") when notebook-local artifacts are required. ₁

The report’s generated outputs include the following tables (a selection, not the full list): ₁

• data/defense_gene_clusters.tsv.gz: 930,573 rows, one per gene-cluster–marker hit, recording system, subtype and source. ₁
• data/species_defense_matrix.tsv.gz: 27,626 rows of species × system presence plus covariates. ₁
• data/species_prophage_burden.tsv.gz: 27,702 rows of species × prophage-module presence plus n_prophage_clusters. ₁
• data/syndrome_pairs.tsv: 28 pairwise co-occurrence rows recording OR, null mean, z, empirical p and BH-FDR q. ₁

The prophage-burden row count differs from the 27,690 pangenome species the report cites elsewhere, and the report does not explain the difference. ₁

Prophage classifier

The eggNOG-description prophage classifier produced 4,005,537 gene-cluster hits, classified into 4,228,150 module-hits across 7 modules. Species counts by number of prophage modules present were: 2 modules, 3 species (<0.01 % of the pangenome); 3 modules, 178 (0.6 %); 4 modules, 4,577 (16.5 %); 5 modules, 6,897 (24.9 %); 6 modules, 6,389 (23.1 %); and 7 modules, 9,658 (34.9 %). ₁

Caveats and Limitations

CRISPR-Cas prevalence (96 %) is inflated by permissive EggNOG description matching and should be read as an upper bound. The Pfam Cas1 (PF01867) hit is present in only 55 % of species. The combined EggNOG-or-Pfam presence set used for syndrome analysis prioritizes recall, so CRISPR-Cas syndrome specificity may be reduced even though the significance results are not expected to be invalidated. ₁

The DISARM accessory-enrichment result is unreliable because PF00176, the DrmB SNF2 helicase Pfam, is a widespread housekeeping-helicase family. Its inclusion inflates the DISARM cluster set with non-DISARM SNF2 hits, which biases the DISARM core/auxiliary/singleton distribution toward core. A PADLOC MacSyFinder-style HMM-plus-context refinement is needed; the report states that arms-race and syndrome results involving DISARM are unaffected because they use species-level presence/absence rather than per-cluster classification. ₁

Retron detection uses the broad RVT_1 Pfam PF00078. The stringent call requires at least one other narrow defense system and is therefore a defense-context filter rather than a retron-specificity filter; it reduced candidates from 15,109 to 15,098, dropping 11 species. Results should be interpreted as reverse-transcriptase candidates in defense-syndrome context, not as characterized retron systems. ₁

The eggNOG-description prophage classifier deliberately captures broad module candidates including integrase, holin, endolysin, CI-like repressor, and tail proteins, which can also match non-phage bacterial genes. Consequently, n_prophage_modules saturates at 7 for 35% of species; the primary arms-race predictor is instead n_prophage_clusters, an unbounded continuous count of matched clusters. n_prophage_modules serves only as a coarse categorical secondary predictor. ₁

The ≥5-genome filter reduced the analysis set to 7,323 of 27,690 pangenome species, or 26%, biasing the analysis toward well-sampled, culturable, high-priority organisms. Extrapolation to the full 293K-genome tree, including many environmental MAGs, remains untested. ₁

The arms-race analysis controls for phylum categorically but does not use a phylogenetically corrected regression. Its partial correlations and negative-binomial GLM are therefore consistent with, but do not formally establish, a phylogenetically independent arms race. ₁

The negative-binomial GLM fixes dispersion at alpha = 1.0 rather than estimating it. Both focal coefficients remained significant at p < 0.001 with wide margins, but standard errors may be slightly under- or over-stated; refitting with an estimated alpha would sharpen inference. ₁

Future Directions

The report recommends PADLOC MacSyFinder-style multi-PFam and gene-order rules for Retron, DISARM and Gabija. It also recommends a GTDB-tree phylogenetic mixed-effects arms-race model. ₁

The report proposes extending the current focused 7-system panel to 20+ defense systems using DefenseFinder/PADLOC HMMs. The proposed additions are CRISPR-Cas subtypes, Zorya, Thoeris, Wadjet, Druantia, pAgo, PARIS, ThsA-ThsB and defense-associated antitoxin cassettes. A larger panel would enable a full defense-repertoire principal-component analysis and genus-level clustering of “defensive phenotype”. This extension is proposed, not completed. ₁

For the R-M Type II × Gabija association (OR = 24), the report proposes mechanistic follow-up. It would first identify species carrying only R-M Type II versus those carrying both systems. It would then run phage-challenge experiments in a tractable host (e.g., E. coli, B. subtilis) to test whether adding Gabija gives a quantitative fitness advantage. That advantage is not established in the report. ₁

Habitat specificity of defense syndromes also remains untested. The report proposes cross-referencing defense repertoires with NCBI-env isolation sources or AlphaEarth embeddings. This environmental layer is a third pillar that the report deferred at the plan-review checkpoint. ₁

Slots Into

• cofitness-network-architecture — the 27-of-28 defense-system co-occurrence pattern and the R-M Type II × Gabija syndrome extend network-level analysis of coordinated gene-system architecture. ₁
• gene-function-acquisition-depth — accessory and singleton enrichment of six defense systems, together with the DISARM detection caveat, provides evidence about mobile acquisition and annotation depth. ₁
• environmental-resistome — the proposed environmental-defense layer and the report’s pan-bacterial defense repertoire provide a basis for testing habitat-specific defense syndromes. ₁
• phage-defense-syndromes-and-arms-race — this concept draws on several of the report’s findings. They include the 7.2%–96.1% prevalence spectrum, the 27-of-28 significant co-occurring pairs led by R-M Type II × Gabija, and the CRISPR-Cas × CBASS null. They also include the phylum-variable prophage–defense association and its framing against Vassallo et al. ₁
• two-speed-bacterial-genome — six of seven defense-system families are enriched in the auxiliary and singleton pangenome relative to the 46.8% core background. ₁
• horizontal-gene-transfer-driven-innovation — elevated singleton fractions are consistent with defense-island horizontal transfer but do not directly demonstrate it. ₁
• scale-dependent-mobile-element-associations — the report infers mobility of defense loci from pan-bacterial singleton enrichment rather than from element-level evidence. ₁
• functional-marker-validation — CRISPR-Cas prevalence depends on the detection method (eggNOG 96% versus ~55% for Cas1 PF01867), and the PF00176 anchor creates a DISARM artefact. The report also contrasts broad-marker with specific-marker calls, and its broad prophage classifier saturates at 7 modules. ₁
• cross-tenant-data-bridging — interproscan_domains and bakta_pfam_domains differ in marker coverage and in their versioned or version-free Pfam accession formats. ₁
• homology-search-negative-evidence — Cas1 PF01867 returns 0 hits in bakta_pfam_domains versus ~25K in InterProScan, a table-coverage false negative. ₁
• sampling-depth-and-downsampling-effects — the ≥5-genome requirement leaves 7,323 species (26 % of 27,690), biasing results toward well-sampled organisms. Extrapolation to the full 293K-genome tree is untested. ₁
• phylogenetic-confounding-of-pangenome-associations — the arms-race models control phylum only categorically, so they are consistent with, but do not formally establish, a phylogenetically independent arms race. ₁