To investigate the diversity and distribution of defense systems, we performed a comparative genomic analysis using complete genome sequences of clinically relevant bacterial species associated with sepsis and other severe infections. Genomes were retrieved from the NCBI GenBank database (as of 2025-12-29) and filtered to include only high-quality, closed, and complete genomes (or chromosomes) to ensure accurate identification and annotation of defense systems. A total of 73 Providencia strains were analyzed: 31 isolates of P. stuartii and 42 isolates of P. rettgeri. All accession numbers for the analyzed genomes are listed in Supplementary Table 1, enabling reproducibility and further validation of our findings. For comparative purposes, we also included complete genomes from 13 additional sepsis-associated bacterial species across diverse taxonomic groups, including Clostridioides difficile, Staphylococcus aureus, Yersinia pestis, Mycobacterium tuberculosis, and Brucella spp., among others. This curated dataset allowed us to systematically compare the repertoire of defense systems across closely related and distantly related pathogens, with a particular focus on the exceptional immune complexity observed in Providencia species.

Statistical comparisons of defense system count between groups were performed using one-way ANOVA with Tukey’s post hoc test (full details in Supplementary Data 1). As illustrated in Figure 1A, Providencia exhibits one of the highest median counts of distinct defense system types among all sepsis-associated bacterial genera analyzed, with a median of approximately 7 systems per genome, exceeding those observed in Brucella, Mycobacterium, and Ralstonia, which show median values below 2 (p < 0.001; see Supplementary Data 1 for full pairwise comparisons). Notably, the distribution of defense system richness in Providencia is not only elevated but also highly variable, as evidenced by the broad interquartile range (IQR) and the presence of several outliers reaching up to 12 distinct systems. In contrast, other clinically relevant genera such as Staphylococcus, Clostridium, and Yersinia display more moderate and less variable defense system profiles, with medians ranging from 3 to 6. At the species level (Figure 1B), this trend is further reinforced: both P. stuartii and P. rettgeri consistently harbor a greater number of defense system types compared to most other sepsis-causing species. Remarkably, while P. rettgeri shows a slightly higher median count than P. stuartii, the two species are nearly indistinguishable in the total number of defense system per genome (p = 0.07), with both exhibiting substantial intra-species variation and frequent occurrences of genomes carrying more than 8 different defense systems. Statistical comparisons reveal significant differences between P. stuartii and P. rettgeri and most other species (p < 0.05), with black asterisks indicating significance relative to P. stuartii and red asterisks indicating significance relative to P. rettgeri.

To further resolve the architectural complexity of these immune systems, we analyzed the abundance of defense system subtypes, specific molecular variants within each system class (e.g., CRISPR-Cas type I-E, II-A; R-M Type I, II, etc.), across the same set of genomes. As shown in Supplementary Figure 1, Providencia again stands out with the highest median subtype count (~8 subtypes per genome), surpassing all other genera (p < 0.001). At the species level (Supplementary Figure 1B), P. stuartii and P. rettgeri exhibit higher subtype richness compared to most other clinical pathogens, with some strains harboring over 10 distinct subtypes. Significant differences were observed between P. stuartii (or P. rettgeri) and most other species (p < 0.05), marked by black and red asterisks, respectively; full statistical details are available in Supplementary Data 1. This elevated subtype diversity indicates not only a greater number of defense systems but also a broader functional repertoire, suggesting that Providencia may deploy a multi-layered, modular defense strategy against invading genetic elements.

To gain deeper insight into the composition and prevalence of defense systems in Providencia, we performed a detailed profiling of immune mechanisms across all analyzed genomes of P. stuartii (n = 31) and P. rettgeri (n = 42). As shown in Figures 2A, B, the most abundant defense system type in both species is the restriction-modification (RM) system, with 61 occurrences (20.75%) in P. stuartii and 77 occurrences (21.04%) in P. rettgeri, respectively, highlighting its fundamental role in innate immunity against foreign DNA. The CRISPR-Cas system is among the most abundant defense systems in both species, detected in 30 strains (10.20%) in P. stuartii and 14 strains (3.83%) in P. rettgeri, representing 10.73% and 11.99% of total defense systems, respectively. Other systems include GAPS2 (30, 10.20% in P. stuartii; 39, 10.66% in P. rettgeri) and Gabija (9, 3.06% in P. stuartii; 10, 2.73% in P. rettgeri), underscoring the presence of multiple non-CRISPR defense mechanisms. Notably, while RM and CRISPR-Cas systems are prevalent in both species, P. rettgeri exhibits higher counts of several non-CRISPR systems such as MazEF (39, 10.66%), Retron (15, 4.10%), and Lamassu-Fam (11, 3.01%), suggesting potential differences in evolutionary adaptation or niche-specific immune pressures. In contrast, P. stuartii shows greater representation of PsyrTA, Mokosh (18, 6.12%), and AbiE (9, 3.06%), indicating distinct immune prioritization between the two species. The full spectrum of defense systems is shown in the main panels (Figures 2A, B), with inset plots detailing the composition of the “Others” category, which accounts for 13.27% in P. stuartii and 11.20% in P. rettgeri. This comprehensive view confirms the presence of numerous rare but functionally distinct systems, including ShosTA, Thoeris, pAgo, and Kiwa, with some occurring in only one or two strains, emphasizing the high degree of genomic heterogeneity and modular evolution within the genus. All raw counts and system classifications are provided in Supplementary Data 2. Together, these findings illustrate a complex, multi-layered defense architecture in Providencia, dominated by RM and CRISPR-Cas systems but enriched by a wide array of auxiliary immune mechanisms that likely contribute to phage resistance and genome stability in diverse environments.

To expand the scope of our defense system profiling beyond complete genomes, we analyzed P. stuartii and P. rettgeri using additional genomic assemblies at the contig and scaffold levels, which are more widely available in public databases like NCBI. Specifically, we applied DefenseFinder to 429 contigs and 50 scaffolds of P. stuartii, and 334 contigs and 244 scaffolds of P. rettgeri. As shown in Supplementary Figures 2A–D, this broader dataset reveals a significantly higher total number of defense system occurrences due to increased genomic sampling. In P. stuartii, the most abundant systems across all contigs are RM (608, 15.31%) and Cas (429, 10.80%), followed by GAPS2 (425, 10.70%) and PsyrTA (409, 10.30%). Similarly, in P. rettgeri, RM (567, 18.50%) and GAPS2 (326, 10.64%) dominate, with MazEF (322, 10.51%) and Gabija (128, 4.18%) also highly represented. Notably, the relative abundance of major systems such as RM, Cas, GAPS2, and PsyrTA remains consistent across complete, contig and scaffold datasets, suggesting robust detection of prevalent defense mechanisms even in fragmented assemblies. While full-length genome assemblies are ideal for accurate defense system annotation, since correct gene order, operon structure, and strand orientation are critical for functional inference, the results from contig and scaffold data demonstrate that they can still capture dominant immune features. This suggests that such fragmented datasets may serve as useful proxies for preliminary assessments of defense system prevalence, particularly when complete genomes are scarce.

We next analyzed the distribution of defense system subtypes, offering a higher-resolution view of functional diversity. As shown in Figures 3A, B, the most abundant subtype across both species is RM subtypes, exhibit substantial diversity, with RM_Type_I, RM_Type_II, RM_Type_IV, and et al. collectively accounting for more than 20% of all subtypes, reflecting a multi-layered DNA surveillance mechanism. Notably, the CRISPR-Cas system is represented almost exclusively by a single subtype: CAS_Class1-Subtype-I-F, across both species, with 30 occurrences (10.22%) in P. stuartii and 14 (3.83%) in P. rettgeri. This subtype dominates the CRISPR-Cas landscape, and no other known CRISPR subtypes such as I-E, II-A, or III-B, were detected in either species, indicating a highly focused evolutionary strategy for adaptive immunity in Providencia. Notably, the “Others” category accounts for 26 (8.84%) in P. stuartii and 54 (14.75%) in P. rettgeri, encompassing rare but functionally distinct variants such as PARIS_I, Thoeris_II, and DS-1, many of which appear in only one or two isolates. The full list of subtypes and their counts is provided in Supplementary Data 2. Together, these findings highlight a unique immune profile in Providencia: dominated by a narrow set of highly prevalent subtypes, yet enriched with a diverse array of auxiliary systems.

To assess the consistency and comprehensiveness of subtype annotation, we performed an independent analysis using the PADLOC tool on the same set of complete genomes, complementing our prior DefenseFinder results. As shown in Supplementary Figures 3A, B, PADLOC identifies a broadly similar profile of dominant defense systems, with RM types collectively being the most abundant subtype in both species. And cas_type_I-F1 also highly represented, confirming the robust detection of these major immune modules across platforms. This concordance supports the reliability of the core defense repertoire observed in our initial analysis. However, PADLOC also reveals several previously undetected or underrepresented systems, particularly within the DMS_other and PDC family, which are not annotated by DefenseFinder. While both tools agree on the dominance of RM and CRISPR-Cas systems, PADLOC’s sensitivity to divergent or atypical architectures allows it to detect additional immune variants, highlighting the importance of multi-tool validation in defense system profiling. However, it should be noted that the PADLOC database has not been updated since its last release, and thus lacks recently characterized systems. Unlike DefenseFinder, which directly annotates both system types and subtypes, PADLOC outputs predictions exclusively at the subtype level, necessitating post grouping to reconstruct higher-order system categories.

To examine the inter-strain variation in immune defense systems, we analyzed the genomic composition and organization of defense system types across a representative subset of P. stuartii and P. rettgeri isolates. Figures 4A, B show stacked bar charts depicting the number and types of defense system subtypes present in individual strains, revealing substantial heterogeneity in immune repertoire composition. For example, P. stuartii strain GCA_010669105 harbors 15 distinct subtypes, while other strains contain as few as 5, indicating significant variability in immune complexity even within the same species. Similarly, P. rettgeri strains exhibit diverse profiles, with some carrying up to 14 subtypes (e.g., GCA_0103188815), while others possess fewer than 4. Stacked bar charts illustrating all analyzed strains can be found in Supplementary Figure 4A (P. stuartii) and 4B (P. rettgeri). The circular genome maps in Figures 4C, D reveal the distribution of defense systems in the genomes of two strains. In P. stuartii GCA_010669105, three major defense-associated genomic islands (GIs) are evident: one near 1.2 Mb; another at ~1.7 Mb; and a third near 3.0 Mb. These clusters are characterized by co-directional gene transcription (indicated by red/blue arrows), consistent with operon-like organization. In contrast, most defense systems appear as isolated, singleton loci in P. rettgeri GCA_010318885; seems to be more randomly located in the genome. The presence of structured defense islands and dispersed systems may suggest a dual evolutionary strategy: one favoring stable, co-adapted gene clusters for robust immunity, and another allowing flexible, piecemeal adaptation to novel threats.

The CRISPR-Cas system is a key component of adaptive immunity in bacteria, and its distribution across Providencia species reveals a striking pattern of both conservation and divergence. As shown in Figures 5A, B, the presence of CRISPR-Cas systems varies significantly between P. stuartii and P. rettgeri. In P. stuartii, the system is nearly ubiquitous, with 30 out of 31 strains (96.77%) harboring at least one CRISPR array, whereas in P. rettgeri, only 15 out of 42 strains (35.71%) possess this defense mechanism, indicating a substantial difference in the prevalence of CRISPR-based immunity between the two species. Actually, this contrast is supported by the subtype-level analysis in Figure 3 and Supplementary Figures 3A, B, which shows that CAS_Class1-Subtype-I-F is the sole CRISPR-Cas subtype detected across all analyzed genomes (All data are provided in Supplementary Data 3). The exclusive dominance of this single subtype suggests a highly specialized evolutionary trajectory for adaptive immunity in Providencia. Notably, while both species carry the same subtype, its relative abundance differs markedly, indicating not only differential acquisition but also varying degrees of functional integration into the immune repertoire. These findings collectively highlight a unique CRISPR landscape in Providencia: characterized by a narrow specificity of Cas subtype and a pronounced species-level disparity in prevalence, suggesting distinct evolutionary strategies for phage defense.

At the genomic level, the architecture of the CRISPR system further distinguishes P. stuartii from P. rettgeri, despite both species typically encoding only a single copy of the core Cas operon. As illustrated in Figures 5C, D, the Cas protein cluster, comprising cas1, cas3, and the other cas genes, is highly conserved in gene content and order across both species, consistent with the canonical Class 1 Type I-F system. However, the organization of CRISPR arrays relative to this Cas locus is different. In P. stuartii, multiple CRISPR arrays are often found dispersed across the chromosome, with some located near the Cas operon and others situated at distant genomic loci, up to several hundred kilobases away. This arrangement suggests that a single Cas complex may function with multiple spatially separated arrays,. In contrast, P. rettgeri strains exhibit a much more restricted configuration: CRISPR arrays are almost exclusively located immediately upstream or downstream of the Cas operon, forming a compact, self-contained unit. This localized architecture implies a more limited and tightly coupled system, where the Cas machinery acts primarily on a single or few adjacent arrays. The distinct genomic distribution of CRISPR repeats thus reinforces the functional divergence between the two species, not only is the CAS_Class1-I-F system more prevalent in P. stuartii, but its genomic integration also supports a more expansive and potentially flexible immune strategy.

The number and distribution of CRISPR arrays not only reflect the complexity of the immune system but may also influence spacer acquisition efficiency and system stability. While the structure of typical CRISPR arrays has been illustrated, whether all CRISPR arrays follow the same trend in terms of their distribution and quantity. To further evaluate the architectural diversity, we analyzed the number of CRISPR arrays per strain across all CRISPR-positive isolates (Figure 5E). The boxplot reveals a difference between P. stuartii and P. rettgeri: while P. stuartii strains harbor 2 to 6 CRISPR arrays, with a median of 5 arrays per genome, P. rettgeri strains are predominantly limited to 2 arrays, with only a few exceptions reaching up to 3. Importantly, this multi-array architecture aligns with the dispersed genomic organization seen in Figure 5C, where multiple CRISPR repeats are located at distant chromosomal positions. In contrast, the limited number of arrays in P. rettgeri, often clustered near the Cas operon, reflects a more streamlined and possibly less flexible immune configuration. Together, these results indicate that the expansion of CRISPR arrays is an important determinant of adaptive immunity capacity, and that P. stuartii has evolved a more complex and scalable system compared to P. rettgeri, reinforcing the notion of species-specific evolutionary trajectories in CRISPR-mediated defense.

Given that differences already exist between P. stuartii and P. rettgeri in terms of CRISPR system distribution and genomic structure, do these two species exhibit further differentiation at the functional level, specifically regarding the number and diversity of spacers. To assess the functional potential, we analyzed the number of spacers per strain across all CRISPR-positive isolates (Figure 5F). The median spacer count in P. stuartii is much higher, approximately 35 spacers per strain, with a broad interquartile range and multiple strains exceeding 60 spacers, indicating not only a greater capacity to record past infections but also a more extensive immunological memory. In contrast, P. rettgeri strains exhibit a substantially lower median of ~25 spacers, with most strains clustering between 20 and 30, and only one outlier surpassing 70. Furthermore, the distribution of spacer counts in P. stuartii is more heterogeneous, whereas P. rettgeri displays a narrower and more uniform profile. This disparity suggests that P. stuartii not only maintains CRISPR-Cas systems at a higher frequency (Figure 5A) and with greater genomic flexibility (Figure 5C), but also accumulates spacers at an elevated rate, likely due to more frequent exposure to phages or more efficient spacer acquisition mechanisms.

CRISPR spacers serve as molecular memory of past encounters between bacteria and MGEs, with each spacer derived from a protospacer sequence in an invading phage or plasmid. By matching these spacers to known MGE databases (the BacMGEnet pipeline), we reconstructed spacer-MGE interaction networks for P. stuartii and P. rettgeri. In total, 546 unique spacers from P. stuartii and 363 from P. rettgeri were identified, of which 110 and 61 matched phages or plasmids, respectively. As shown in Figures 6A, B, the resulting networks reveal distinct patterns of spacer-MGE interactions. The green arrows in the P. stuartii spacer-MGE interaction network highlight two distinct clusters where multiple spacers from different strains converge on a single phage (yellow oval), indicating shared immunity against common viral threats. These highly connected hubs suggest that certain phages have repeatedly infected multiple P. stuartii strains, leading to the acquisition of identical or highly similar spacers. In contrast, the P. rettgeri network lacks such densely interconnected MGE nodes, with most phages linked to limited strains. These findings suggest that CRISPR-based immunity plays a key role in shaping bacterial population dynamics.

In addition, spacers enable the reconstruction of host-MGE interaction networks by linking bacterial strains to the phages or plasmids they have encountered. In our analysis, the P. stuartii network comprises 56 nodes: 18 host strains and 38 MGEs (37 phages and 1 plasmid); while the P. rettgeri network includes 40 nodes: 13 hosts and 27 MGEs (24 phages and 3 plasmids) (Supplementary Figures 5A, B). The P. stuartii network displays a highly interconnected architecture with numerous edges, indicating frequent and shared targeting of phages across multiple strains. By contrast, the P. rettgeri network is markedly sparser, with fewer connections and more isolated host-MGE pairs, consistent with its lower CRISPR array abundance and spacer diversity observed earlier. Notably, several P. stuartii strains interact with multiple phages, and certain phages are targeted by multiple hosts, suggesting the emergence of community-level immunity. Together, these distinctions underscore the value of network-based approaches in uncovering ecological and evolutionary dynamics of host-phage interactions and provide a foundation for future investigations into phage resistance mechanisms and the potential for precision phage therapy in Providencia infections.

Furthermore, the sequence characteristics, diversity, and structural stability of CRISPR repeats are not only the foundation of CRISPR system function but may also influence crRNA processing efficiency, Cas protein recognition ability, and overall immune activity. Analysis of repeat characteristics across Providencia genomes revealed diversity within both species (Table 1). A total of 19 distinct repeat variants were identified, indicating that both species maintain multiple repeat types, potentially supporting functional heterogeneity among CRISPR loci. Notably, the number of spacers associated with each repeat type varies dramatically, from as few as 4 to over 300, reflecting differential expansion dynamics and highlighting that certain repeat loci are hotspots for spacer acquisition and immune memory accumulation. The predicted RNA secondary structures of representative repeats (Supplementary Figures 6, 7) further reveal their molecular architecture. Both P. stuartii and P. rettgeri repeats form stable stem-loop structures, consistent with the canonical hairpin required for Cas protein recognition and crRNA maturation.

In addition to immune defense systems, the pathogenic potential of Providencia species is critically shaped by two key genomic determinants: antibiotic ARGs and VFs. ARGs enable these bacteria to survive under antimicrobial pressure, posing significant challenges in clinical settings where treatment options may be limited. Meanwhile, VFs, such as adhesins, toxins, secretion systems, and iron acquisition systems, facilitate host colonization, tissue invasion, and immune evasion, thereby driving infection progression. The co-occurrence of defense systems with ARGs and VFs within the same genomic contexts can promote the coordinated dissemination of multiple adaptive traits through HGT, potentially giving rise to multidrug-resistant, highly virulent clones. Therefore, understanding the prevalence, diversity, and genomic localization of ARGs and VFs is essential for comprehensively characterizing the evolutionary and clinical significance of Providencia pathogens.

To investigate whether defense systems are associated with antimicrobial resistance and virulence, we performed correlation analyses between the total number of defense system types, CRISPR spacer count, ARG count, and VF count in P. stuartii and P. rettgeri, respectively (Figures 7A, 8A, all data are provided in Supplementary Data 4). The results reveal species-specific patterns. In P. stuartii, a positive correlation is observed between the number of defense system types and ARGs (r = 0.59), suggesting that strains with more diverse immune repertoires tend to harbor a greater number of ARGs. The association is further described by a negative correlation between defense system diversity and VF count (r = −0.82), indicating a potential trade-off between immune complexity and virulence gene acquisition. Notably, CRISPR spacer count shows a positive correlation with VF count (r = 0.90). In contrast, P. rettgeri exhibits a weaker but still positive correlation between defense types and ARGs (r = 0. 38), while showing a negative correlation with VFs (r = −0.50). These findings suggest that although both species show some degree of integration between defense and resistance traits, P. stuartii displays a more pronounced pattern of co-enrichment of defense and resistance, coupled with a potential antagonism between immunity and virulence.

In Figures 7B, C, we present the correlation analysis between some defense systems identified in P. stuartii and specific ARGs or VFs. Similarly, Figures 8B, C depict these relationships for P. rettgeri. Each point in the figures represents a pairwise comparison of a particular defense system with an individual ARG or VF, using the Phi coefficient to quantify the strength and direction of association. For P. stuartii, several defense systems exhibit strong correlations with certain ARGs, suggesting potential co-occurrence patterns within the genome. Conversely, a few defense systems display negative correlations hinting at possible antagonistic relationships or genetic linkage constraints. In P. rettgeri, the overall pattern is less pronounced but follows a similar trend, with weaker but still observable correlations between selected defense systems and both ARGs and VFs. These findings provide insight into the genomic architecture of defense, resistance, and virulence traits, highlighting the complex interplay among them.

While the pairwise correlation analyses reveal several statistically significant associations, these findings should be interpreted with caution due to underlying genomic architecture and evolutionary dynamics. Notably, ARGs are frequently observed in clusters, often co-occurring within MGEs such as plasmids, transposons, or integrative and conjugative elements. Consequently, the observed correlations between defense systems and ARGs could reflect genomic co-localization on shared mobile platforms, rather than a direct functional or selective linkage. Similarly, the negative correlation between defense systems and VFs may arise from genomic space constraints, fitness trade-offs, or differential niche adaptation, rather than direct antagonism. Therefore, while the statistical associations provide valuable hypotheses about potential interactions among defense, resistance, and virulence traits, they likely capture indirect signals shaped by genome plasticity and mobile element dynamics. To determine whether these correlations have true biological significance, future studies should integrate genomic context analysis, transcriptomic data, and experimental validation.

The evolution of defense systems in bacteria is a dynamic process marked by both the gain and loss of gene families, which are critical for survival against foreign genetic elements such as bacteriophages and plasmids. Previous studies have highlighted that the acquisition of new defense mechanisms often occurs through HGT, enabling rapid adaptation to environmental pressures. Conversely, the loss of certain defense-related genes may be driven by fitness costs associated with maintaining these systems or due to changes in selective pressures over time. In this study, we calculated the gain and loss events within defense system gene families by employing comparative genomic analyses and ancestral state reconstruction methods. By understanding the dynamics of gene gain and loss, we can better comprehend how bacteria evolve their defensive capabilities in response to changing environments, providing insights into the complex interplay between microbial communities and their biotic and abiotic surroundings.

To investigate the evolutionary dynamics of defense systems in Providencia, we reconstructed gene gain and loss events across the phylogeny of P. stuartii and P. rettgeri using a phylogenomic framework. OrthoFinder was first employed to cluster protein-coding genes into orthogroups, enabling the inference of gene family presence/absence profiles across strains. Subsequently, the COUNT software was applied under a maximum-likelihood birth-death model to estimate per-gene-family rates of gain and loss, accounting for lineage-specific evolutionary processes. This approach allows for the modeling of gene content evolution by treating gains and losses along the branches of the species tree (All data are provided in Supplementary Data 5).

Our results reveal distinct evolutionary trajectories between defense systems and the whole genome in both species. In P. stuartii, defense systems exhibit experience a higher rate of gene gain (mean: 1.90 × 10^-2) compared to the whole genome (1.10 × 10^-2), yet they also undergo substantial gene loss (3.78 × 10^-2), which—although slightly lower than the whole-genome loss rate (4.91 × 10^-2) (Table 2), resulting in a net loss trend. This suggests that while P. stuartii actively acquires new defense genes, it simultaneously discards others at a high rate, indicative of a dynamic, high-turnover immune arsenal. In contrast, P. rettgeri displays a higher average gain rate in defense systems, driven by an elevated gain rate (5.01 × 10^-2), nearly double its whole-genome gain (2.62 × 10^-2), and a low loss rate (5.96 × 10^-2), which is only about one-third of its genomic loss rate (1.74 × 10^-2) (Table 3). This indicates that, unlike P. stuartii, P. rettgeri is undergoing expansion and stabilization of its defense repertoire, possibly reflecting stronger selective pressures to maintain diverse immunity mechanisms. Collectively, these findings demonstrate that defense systems in Providencia are subject to high evolutionary turnover, but with species-specific strategies: P. stuartii maintains a dynamic, high-turnover, replaceable defense repertoire, whereas P. rettgeri favors net acquisition and retention of defense genes. This divergence highlights the distinct evolutionary paths taken by these closely related pathogens in adapting their immune arsenals.

To experimentally test whether these computationally identified systems are functional, we selected two well characterized systems, Gabija from P. stuartii (GCA_010669105) and Septu from P. rettgeri (GCA_010318885), for heterologous expression and phage challenge assays in E. coli DH5 alpha (Supplementary Data 6).

Gabija and Septu were chosen for validation due to their well-characterized molecular mechanisms. Both systems have been extensively studied in E. coli, where they confer robust immunity against diverse phages through distinct abortive infection or DNA-targeting strategies. As shown in Figure 9A, the wild-type Gabija system from P. stuartii strongly inhibited plaque formation by T4 phage, confirming its anti-phage activity. Introducing a single point mutation (E465K) in the GajA subunit, a residue conserved in the predicted ATPase domain, significantly attenuated this defense. Similarly, the Septu system from P. rettgeri conferred potent resistance against T7 phage (Figure 9B), while the PtuB H53K mutant lost nearly all protective capacity. Broader phage profiling revealed that Gabija exhibited strong activity against T4 (+++), weak inhibition of T7 (+), and no effect on λ phage (−), whereas Septu was highly effective against T7 (+++) but inactive against both T4 and λ (Figure 9C). These results demonstrate that Providencia-derived defense systems retain functionality in a heterologous host and display phage -specific activity patterns consistent with their known biological roles.

Notably, although Providencia belongs to the Enterobacteriaceae family and shares close phylogenetic proximity with E. coli, the Gabija system from P. stuartii GCA_010669105 exhibits structural divergence. While core domains (e.g., TOPRIM in GajA) are conserved, the GajA protein is extended to 736 amino acids, with a longer C-terminal region than the well characterized homolog; however, we do not know whether the C-terminal has any function or not. This uncharacterized segment may reflect Providencia-specific adaptations. However, comprehensive functional dissection of such features requires phages naturally infecting Providencia, which are currently unavailable for us.

All complete genome sequences of P. stuartii and P. rettgeri were retrieved from the NCBI GenBank database (https://www.ncbi.nlm.nih.gov/genbank/) in 2025-12-29. A total of 31 P. stuartii and 42 P. rettgeri isolates with complete genomes or chromosome available were included in this study. To broaden our analysis beyond complete genomes, we retrieved all publicly available genomic assemblies from NCBI, including contig- and scaffold-level sequences. This yielded 429 contigs and 50 scaffolds for P. stuartii, and 334 contigs and 244 scaffolds for P. rettgeri. All assemblies were processed with DefenseFinder under default parameters to identify and classify defense systems. Only complete genomes were selected for most analysis to ensure high genomic integrity and to minimize assembly-related artifacts in the identification and genomic context analysis of defense systems. All accession numbers for the analyzed genomes are provided in Supplementary Table 1. For downstream analyses requiring different genomic formats (e.g., nucleotide sequences in FASTA, annotated features in GFF or GBK), corresponding files were downloaded directly from NCBI without re-annotation to maintain consistency with the reference records and to ensure reproducibility of the analyses. However, due to differences in file availability across NCBI repositories and formats, slight discrepancies in isolate representation may occur between individual bioinformatics workflows.

Defense systems were identified using DefenseFinder (version 2.2) with default parameters, a widely used tool for the detection and classification of known antiviral defense systems in prokaryotic genomes (27, 58). The analysis was performed on the complete genome sequences of 31 P. stuartii and 42 P. rettgeri isolates. DefenseFinder was run for each isolate, and results were aggregated for comparative analysis. To visualize the distribution and composition of defense systems across strains, Python (v3.9) libraries such as pandas, matplotlib, and seaborn are used for data processing and plotting.

A complementary analysis using PADLOC (v2.0.0) on the same set of complete P. stuartii and P. rettgeri genomes previously analyzed with DefenseFinder (59, 60). PADLOC was run under default settings to identify defense systems at the subtype level. Results from PADLOC were compared with DefenseFinder outputs to assess concordance in major system detection and to identify potentially missed or underrepresented immune modules.

The CRISPR-Cas systems were identified using CRISPRCasFinder (v4.2.21) (61) and CRISPRCasTyper (v1.3) (62) with default parameters, run locally in command-line mode to enable batch processing. CRISPRCasFinder was used to detect both CRISPR arrays and associated cas genes across all complete genomes. Only CRISPR arrays with an evidence level of 3 or 4, indicating strong support based on repeat conservation and spacer length, were retained for downstream analysis to ensure high-confidence predictions. Genomic context analysis was performed to determine the physical linkage between CRISPR arrays and cas gene clusters. Manual curation was carried out using genome browsers and sequence inspection to confirm the boundaries of arrays and the orientation of flanking genes. The RNA secondary structures of the repeat sequences were predicted using the RNAfold WebServer (available at http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) with default settings. The minimum free energy (MFE) structure was used to assess the potential for stable stem-loop formation.

To investigate the potential interactions between P. stuartii and P. rettgeri with MGEs (including phages and plasmids), the BacMGEnet computational pipeline (https://github.com/mgtools/BacMGEnet) was employed in this study (63). Briefly, the downloaded genomic data of the two bacterial species was used as input. CRISPR arrays and associated spacers were identified using the integrated method, followed by redundancy removal with CD-HIT-EST to obtain unique spacers. These unique spacers were then queried against the multi-source MGE database (downloaded from BacMGEnet pipeline) integrated by BacMGEnet (encompassing phage databases such as GPD, MVP, RVDB, and mMGE, as well as plasmid databases including COMPASS and PLSDB) using BLASTN. Matches were filtered with criteria of >90% sequence identity, >80% query coverage, and e-value <0.001 to identify potential MGEs containing corresponding protospacers. A greedy algorithm was applied to select non-redundant bacterial hosts and MGEs, which were used to construct Host-MGE and Spacer-MGE interaction networks. Network visualization along with manual verification were conducted using Cytoscape.

To investigate the potential associations between defense systems and genomic features related to pathogenicity, VFs and ARGs were annotated for each Providencia genome using data retrieved from the gcPathogen database (https://nmdc.cn/gcpathogen/), a curated resource for pathogen-associated genes in bacterial genomes (64). The presence and abundance of VFs and VRGs were matched to the corresponding isolates based on genome accession ids.

A pairwise correlation analysis was performed between the occurrence profiles of defense system and the presence of individual VFs and ARGs. The analysis was conducted using Phi coefficient (φ), a measure of association for binary variables. To evaluate the association between the overall abundance of defense systems and ARGs/VFs, Spearman’s rank correlation coefficient (ρ) was calculated to assess non-linear relationships, and the corresponding p-values were adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR). Correlations with an adjusted p-value < 0.1 were considered suggestive of potential biological associations.

To infer the evolutionary dynamics of defense system gene families in Providencia, we applied a probabilistic birth-and-death model implemented in COUNT v10.04 (https://www.iro.umontreal.ca/~csuros/gene_content/count.html) (65). The analysis was based on defense system-specific gene family clustering and a species tree inferred using OrthoFinder (v2.5.4) with default parameters (66). OrthoFinder was used to identify orthogroups from the proteomes of Providencia isolates (31 P. stuartii and 42 P. rettgeri strains available from NCBI GenBank), and the rooted species tree generated by OrthoFinder served as the phylogenetic framework for downstream ancestral state reconstruction. The gene family gain, loss, and duplication rates were estimated using the gain-loss-duplication model under a Poisson distribution, with rate variation across branches modeled using three discrete categories of the gamma distribution to account for heterogeneity in evolutionary rates. Similar results were obtained when using four gamma categories, confirming the robustness of the model. For optimization, 100 independent runs were performed to ensure convergence, with the likelihood threshold set at 0.1. Custom scripts were written in-house to quantify gains and losses of defense system families using COUNT-derived results; these scripts are available in the Supplementary Materials (Supplementary Data 7). This combined OrthoFinder-COUNT approach enabled the reconstruction of ancestral gene content and the estimation of lineage-specific rates of defense system gene family evolution across the Providencia phylogeny.

The Gabija and Septu anti-phage defense systems were identified in representative Providencia strains: Gabija from P. stuartii (GCA_010669105) and Septu from P. rettgeri (GCA_010318885). Full genomic coordinates, gene accession IDs, protein sequences, and complete DNA sequences of both systems are provided in Supplementary Data 6. For functional assays, the entire defense loci, including 500 bp upstream and 200 bp downstream of the coding regions, were chemically synthesized by Sangon Biotech Co., Ltd. (Shanghai, China) and cloned into the low-copy-number plasmid vector pBR322 under native regulatory elements.

Site-directed mutagenesis was performed to generate loss-of-function variants: GajA E465K in the Gabija system and PtuB H53K in the Septu system. Mutations were introduced using the QuickMutation Site-Directed Mutagenesis Kit (Cat. No. D0206S, Beyotime Biotechnology, China), following the manufacturer’s protocol. All constructs were verified by Sanger sequencing.

Wild-type and mutant defense constructs were transformed into Escherichia coli DH5α cells for phenotypic characterization. Phage plaque assays were conducted as follows: overnight cultures of DH5α harboring each construct were diluted 1:100 in LB medium and grown to mid-log phase (OD600 ≈ 0.5). Cells were mixed with serial dilutions (10⁻¹ to 10⁻⁷) of laboratory-maintained bacteriophages T4 (Myoviridae), T7 (Podoviridae), or λ (Siphoviridae), incubated at 37 °C for 15 min to allow phage adsorption, and then drop plated in agar overlays on LB agar plates. Plaques were visualized after overnight incubation at 37 °C. Defense activity was presented based on plaque formation relative to empty-vector controls: “−” indicates no restriction (plaque morphology identical to control), “+” denotes moderate reduction in plaque number or size, and “+++” represents strong inhibition (few or no plaques even at low phage dilutions).

Statistical comparisons of defense system type/subtype counts across bacterial genera and species were conducted using one-way analysis of variance (ANOVA) to assess overall differences between groups. Where significant main effects were detected, Tukey’s honest significant difference (HSD) post hoc test was applied to perform pairwise comparisons between Providencia (genus level), P. stuartii, P. rettgeri (species level), and other target taxa. Significance thresholds were set at p < 0.05 and p < 0.001. All statistical analyses were implemented in Python using the scipy and statsmodels packages. Full pairwise comparison results are provided in Supplementary Data 1.