This study was conducted in accordance with the ethical standards of Hospital Universitario Ramón y Cajal. Fifty anonymized, non-diarrheal stool samples previously submitted to the Microbiology Department for routine parasitological screening were collected without identifiable patient information. Samples were pooled into five groups of ten individual specimens each to facilitate BALO enrichment and detection.

B. bacteriovorus HD100 (ATCC 15356) was used as control predator strain and Pseudomonas putida KT2440 as the reference prey (Martínez et al., 2016). Additional prey strains included P. aeruginosa PAO1 (ATCC 47085) and Escherichia coli (ATCC 25922), and clinical isolates of P. aeruginosa obtained from cystic fibrosis patients (Díez-Aguilar et al., 2021), and Klebsiella pneumoniae isolates from urinary tract infections (Avendaño-Ortiz et al., 2023). Clinical strains were selected based on their antibiotic resistance profiles and previously documented lipopolysaccharide modifications. The strains used in this study are listed in Supplementary Table S1.

The prey strains were cultured from glycerol stock in blood agar (CAN-I-BD) medium at 37 °C and were cultivated in liquid with shaking at 180 rpm in LB medium. All predation co-cultures were performed in diluted nutrient broth (DNB, 1:10 dilution of standard nutrient broth) supplemented with 2 mM CaCl₂ and 1 mM MgSO₄ (Herencias et al., 2017).

BALO enrichment was performed following an adaptation of established protocols (Iebba et al., 2013; Herencias et al., 2017). Ten grams of pooled fecal material were suspended in 50 mL of sterile saline solution (0.9% NaCl) and vortexed vigorously for 5 min to homogenize. To this suspension, 5 mL of DNB containing 10⁸ CFU/mL of P. putida KT2440 was added as prey. The enrichment culture was incubated at 30 °C with shaking at 150 rpm for 48 h to allow the growth of predatory bacteria.

Following incubation, cultures were filtered through 0.45 μm cellulose acetate membrane filters (Millipore) to remove most intact prey cells while allowing small predatory bacteria to pass. One milliliter of filtrate was transferred to tubes containing 500 μL of each potential prey strain adjusted to 10⁸ CFU/mL and co-incubated for 24 h at 30 °C and shaking 150 rpm.

BALO viability and enumeration were determined using the double-layer agar overlay method (Bolyen et al., 2019). Briefly, 1.5 mL of co-culture was added to 4 mL of melted DNB soft agar (0.7% w/v), maintained at 50 °C, supplemented with 2 mM CaCl₂ and 1 mM MgSO₄. The mixture was gently vortexed and immediately poured onto plates containing DNB agar (1.5% w/v). Plates were allowed to solidify at room temperature for 30 min before incubation at 30 °C for 48–96 h. Clear lytic plaques as indicative of predatory activity were counted and expressed as plaque-forming units per milliliter (PFU/mL).

Genomic DNA was extracted from fecal pools and purified lytic plaques using the QIAamp DNA Stool Mini Kit (Qiagen, Germany) following the manufacturer’s protocol. For fecal samples, we used an aliquot of 3 mL, centrifuged and proceeded with the extraction protocol. For lytic plaque samples (Pool 1, 2, 4, 5), plaques were excised from soft agar, resuspended in 500 μL of sterile phosphate-buffered saline, and centrifuged at 10,000 × g for 10 min to pellet cells before proceeding with the extraction protocol. DNA concentration and purity were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific), and DNA integrity was evaluated by agarose gel electrophoresis. Extracted DNA was stored at −20 °C until further analysis.

To confirm the presence of predatory bacteria, PCR amplification targeting the Bdellovibrionaceae family was performed using specific primers Bd347F (5′-ATAAGGGATGACGACGACGGAGG-3′) and Bd549R (5′-GCTAG GATCCCTCGTCTTACC-3′), following the conditions previously established (Iebba et al., 2013).

To characterize the bacterial composition of fecal pools and lytic plaques, 16S rDNA gene amplicon sequencing targeting the V3-V4 hypervariable region was performed. PCR amplification was conducted using universal primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′) with Illumina adapter sequences. Amplicons were sequenced on an Illumina MiSeq platform using the MiSeq Reagent Kit v3 (2 × 300 bp) for 600 cycles.

Raw sequencing reads were processed using the QIIME 2 bioinformatics platform (Bolyen et al., 2019). Quality filtering, denoising, and ASV identification were performed using the DADA2 (Callahan et al., 2016). Forward and reverse reads were trimmed to remove primers and low-quality bases, and reads were merged with a minimum overlap of 12 bp. Chimeric sequences were identified and removed using the consensus method in DADA2. Taxonomic classification of ASVs was performed using the SILVA 138 database (Quast et al., 2013) with the classify-sklearn naive Bayes classifier. ASVs with total abundances below 10 reads across all samples were removed to control for spurious sequences and sequencing errors.

Due to the low abundance of predatory bacteria in co-cultures with prey, DNA enrichment was necessary prior to whole-genome sequencing. Isolated BALOs from pools 1 and 2 were subjected to sequential enrichment by co-culturing with increasing volumes of P. putida prey. Specifically, predators were incubated with 1 mL of prey (10⁸ CFU/mL) in 8 mL HEPES buffer (25 mM, pH 7.8) at 30 °C with shaking at 170 rpm for 72 h. Cultures were centrifuged at 15,000 × g for 5 min at room temperature, supernatants were discarded, and pellets were resuspended in double the previous prey volume with fresh HEPES buffer. This process was repeated four times until sufficient biomass was obtained, as evidenced by visible pellet formation.

Total genomic DNA was extracted using the QIAamp DNA Mini Kit (Qiagen) and quantified using the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). For short-read sequencing, DNA libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina) and sequenced on an Illumina MiSeq platform with the MiSeq Reagent Kit v2 Nano (2 × 150 bp) for 300 cycles, generating 250,000–500,000 reads per sample. Genomic DNA were sequenced with Oxford Nanopore technology at Plasmidsaurus (Eugene, United States¹) using the V10 chemistry library prep kit with the R10.4.1 flow cells. Sequences were annotated using BAKTA (Schwengers et al., 2021).

The 162 deposited genomes were initially annotated using Bakta 1.7.0 (Schwengers et al., 2021). Subsequently, core and accessory genomes were defined using the PATO tool (Fernández-de-Bobadilla et al., 2021). This allowed for an analysis of the population structure, the annotation of adaptive features, and the creation of gene networks. Specifically, we utilized the PATO functions “core_plots” to determine the size of the pangenome, core, and accessory genome, and “similarity_tree” to generate pseudo-phylogenetic trees.

We assigned the species to each isolate using the function “classifier” of PATO (Fernández-de-Bobadilla et al., 2021). “Classifier” assigns each query genome to the closest reference genome from the NCBI² by calculating the ANI of each genome to the reference ones. It assigns a reference species if the ANI is over 95% of identity. The heatmap from Supplementary Figure S3 that included 162 sequenced genomes of Bdellovibrio from the NCBI was created using the function “similarity_tree” and the values from the function “mash” from PATO. Clades were defined based on minimum within-node pairwise ANI scores using the unsupervised clustering MClust (Fraley and Raftery, 2007).

Predatory capacity (PFU counts) data were analyzed using a negative binomial generalized linear model (GLM) to account for overdispersion and zero-inflation inherent in count data. The model included the predator strain and prey type as fixed effects. Pairwise comparisons between isolates (BD_H1, BD_H2) and the reference strain (HD100) were performed using Tukey’s Honestly Significant Difference (HSD) post-hoc tests. Statistical significance was set at p < 0.05. All analyses were conducted using R software (v.4.3.1).

Assembled genome sequences for isolates BD_H1 and BD_H2 have been submitted to the EBI database under the project accession number PRJEB102169.

Four out of the five pooled fecal samples analyzed (pools 1, 2, 4, and 5) yielded lytic plaques following enrichment with P. putida KT2440 as prey bacteria. These plaques displayed characteristic features of BALO predation: round morphology with sharp boundaries, progressive size increase over 3–5 days, and absence of bacterial colonies at their centers (Figure 1A; Supplementary Figure S1A). Viable predatory bacteria were quantified using the double-layer agar overlay method.

16S rDNA gene amplicon sequencing of DNA extracted from lytic plaques revealed variable relative abundances of Bdellovibrionaceae across plaque-derived samples: 28% in pool 1, 40% in pool 2, 5% in pool 4, 38% in pool 5, and 50% in the B. bacteriovorus HD100 positive control (Figure 1B; Supplementary Figure S1B). The remaining bacterial composition consisted primarily of prey organisms (Pseudomonas spp.) and commensal bacteria from the fecal matrix. Notably, despite comprehensive 16S rDNA gene sequencing of the five original fecal pools prior to enrichment, no Bdellovibrionaceae sequences were detected above the detection threshold (>10 reads per amplicon sequence variant, ASV), indicating extremely low natural abundance in the gut microbiota (Supplementary Figure S2). This finding is consistent with previous molecular detection studies showing transient or rare presence of BALOs in human-associated environments (Iebba et al., 2013; de Dios Caballero et al., 2017). Finally, family-specific PCR amplification using primers Bd347F/Bd549R confirmed the presence of Bdellovibrionaceae DNA in Pool1 and Pool2 plaque-positive pools, yielding the expected 202 bp amplicon (Iebba et al., 2013).

After successful enrichment in liquid medium and repeated and robust formation of lytic plaques, two of the predator isolates were selected and sequenced (from pool 1 and pool 2, hereafter BD_H1 and BD_H2, respectively). To establish taxonomic identity and phylogenetic relationships, we employed a comprehensive comparative genomic approach using the PATO pipeline (Fernández-de-Bobadilla et al., 2021), which integrates multiple analytical modules for bacterial genome characterization. Taxonomic assignment was performed using the PATO “classifier” function, which determines species identity by calculating average nucleotide identity (ANI) values against reference genomes in the NCBI database. The classifier assigns a reference species when ANI exceeds 95% identity. Both isolates were confirmed as B. bacteriovorus, with ANI values of 99.10% (BD_H1) and 99.09% (BD_H2) relative to the reference strain HD100 (Supplementary Table S2). These high ANI values, well above the 95% species delineation threshold, confirm conspecific status while indicating strain-level genetic variation.

The BD_H1 genome was assembled into a single circular chromosome of 3,782,514 bp with a GC content of 50.65%. The BD_H2 genome comprised three contigs, identified through hybrid assembly using Illumina and Oxford Nanopore sequencing technologies. The main chromosome (Contig_2, 3,782,519 bp) has a GC content of 50.67%, consistent with BD_H1 and typical of B. bacteriovorus. The additional contigs (Contig_1, 5,354 bp; Contig_3, 5,667 bp) harbor ribosomal RNA operons and transfer RNA genes with elevated GC content (overall assembly: 56.69%), likely representing duplicated chromosomal regions or horizontally transferred elements rather than separate replicons. Detailed analysis of the smaller contigs revealed important structural and functional features. Both Contig_1 and Contig_3 harbor ribosomal RNA (rDNA) operons and transfer RNA (tRNA) genes (Supplementary Table S3), which are typically chromosomally encoded in bacteria.

Comparative alignment showed that Contig_1 shares 100% nucleotide identity with a specific region of the main chromosome (Contig_2) (Supplementary Table S3), strongly suggesting it represents a duplicated chromosomal segment or assembly artifact arising from repetitive rDNA operons rather than an independent replicon. In contrast, Contig_3 exhibited 99.75% nucleotide identity with Stenotrophomonas maltophilia genomic sequences (Supplementary Table S3). The presence of complete rDNA and tRNA genes on this contig, combined with its foreign origin, suggests either: (i) acquisition of a mobile genetic element carrying accessory metabolic functions during gut colonization or (ii) a genuine integrative element that has mobilized between distantly related species. Despite the presence of these additional contigs, the core genome of BD_H2 (Contig_2) exhibited typical B. bacteriovorus architecture and completeness metrics.

To position our human-derived isolates within the broader diversity of predatory bacteria, we retrieved 162 publicly available genomes annotated as Bdellovibrio spp. from the NCBI RefSeq and GenBank databases (Supplementary Table S4). Given the limited number of Bdellovibrio genomes historically available in public databases and the recent taxonomic reclassification of predatory bacteria within the phylum Bdellovibrionota (Waite et al., 2020; Davis et al., 2024), we performed comprehensive genome-wide validation to confirm that deposited sequences truly represent authentic members of the genus B. bacteriovorus. To accomplish this, we employed the PATO “classifier” function, which establishes species boundaries using average nucleotide identity (ANI) calculations against rigorously validated reference genomes.

Pairwise ANI calculations between Bdellovibrio genomes (including our isolates BD_H1, BD_H2) revealed striking phylogenetic heterogeneity within publicly available sequences (Supplementary Figure S3). Most of the deposited genomes showed low ANI values (<90%), indicating they represent distinct species or even genera within the Bdellovibrionaceae family. This finding highlights the substantial taxonomic diversity that remains poorly characterized within this predatory bacterial group and underscores the need for systematic genomic-based reclassification of historical isolates.

Applying the widely accepted >94% ANI threshold for bacterial species delineation (Jain et al., 2018), we identified only 16 genomes (10.4% of total) that qualified as B. bacteriovorus sensu stricto. To examine fine-scale genomic variation within this species complex, we performed a single nucleotide polymorphism (SNP) analysis across the core genome and calculated SNP density normalized per megabase (SNPs/Mb) for all pairwise comparisons (Figure 2A; Supplementary Table S5). SNP-based heatmap analysis revealed distinct patterns of genomic divergence within the B. bacteriovorus clade. BD_H1 and BD_H2 exhibited remarkably low SNP densities when compared to each other and to the reference strain HD100 (Supplementary Table S2). These extremely low SNP densities indicate a recent common ancestry and minimal evolutionary divergence between human-derived isolates and well-characterized environmental strains.

To further resolve the relationships within the B. bacteriovorus sensu stricto and its closest relatives, we selected the 16 genomes with ANI ≥ 95% and an additional 24 genomes with ANI of 90–94%, for a total of 22 high-quality genomes (Supplementary Figure S3). We then performed a core genome analysis and constructed a matrix of pairwise average nucleotide identity (ANI) values, alongside single nucleotide polymorphism (SNP)-based phylogenetic analysis of the core annotated genome (Supplementary Figures S4, S5).

Complementing the ANI-based clustering, the SNP-based maximum likelihood phylogenetic tree of the core genome (Figure 2B) provides additional resolution of relationships among these 22 genomes. The human-derived isolates clustered together with HD100 and closely related reference strains, supporting both their classification and their close evolutionary relationship, while other clades showed longer branches, indicative of additional divergence.

We then analyzed the pangenome composition and evolution across the 22 Bdellovibrio genomes, stratified by core genome, accessory genome, and total pangenomes (Figure 3). The core genome (green line), representing genes present in ≥95% of the genomes, remained remarkably stable and nearly constant across all 22 analyzed strains, plateauing at approximately 2,500–2,650 orthologous gene families. This genomic invariance, even across geographically and ecologically diverse isolates (including human-derived, environmental aquatic, soil, and wastewater strains), indicates strong functional constraints on the predatory lifestyle. The obligate predatory strategy, which requires specialized systems for prey sensing, attachment, invasion, and intracellular lysis, appears to impose a well-defined minimal genomic blueprint that cannot be substantially altered without compromising fitness. In contrast, the accessory genome (red line) exhibited substantial variation and consistent expansion as more genomes were sequenced, reaching approximately 4,100–4,200 gene families by the 22 genomes. The steep trajectory of the red curve indicates that each newly added genome contributes an average of ~60–100 novel genes not found in the previously sequenced strains. Finally, the total pangenome (blue line) demonstrated logarithmic growth, increasing from ~4,000 genes in the first 10 genomes to ~6,500–7,000 genes at n = 22, approaching saturation but not yet fully converging. This incomplete saturation suggests that additional Bdellovibrio diversity remains undescribed, consistent with our finding that only 10.4% of deposited Bdellovibrio-labeled genomes are truly B. bacteriovorus sensu stricto. The continued accumulation of novel accessory genes underscores the substantial hidden diversity within this genus and supports the need for systematic taxonomy-guided exploration and characterization of predatory bacterial isolates from undersampled ecological niches.

The remaining 122 genomes (74.8%) displayed ANI values below 90%, confirming their assignment to divergent lineages within the broader Bdellovibrionaceae family. This distribution reveals that current NCBI annotations significantly oversimplify the taxonomic complexity of predatory bacteria, with many genomes labeled as “Bdellovibrio sp.” actually representing phylogenetically distant taxa. Our analysis demonstrated that the genus Bdellovibrio encompasses far greater genomic diversity than previously recognized, with multiple cryptic species awaiting formal description.

To assess prey range and specificity, both human-derived BD_H1 and BD_H2 predators and the reference strain B. bacteriovorus HD100 were challenged against a panel of clinically relevant prey organisms in three to five independent replicate predation assays (Table 1; Supplementary Table S6).

Predation assays revealed marked strain-specific prey preferences among the three predators tested. The reference strain HD100 demonstrated the broadest prey spectrum, efficiently preying on P. aeruginosa, P. putida, E. coli, and select K. pneumoniae isolates. In contrast, BD_H1 and BD_H2 predators exhibited narrower prey ranges, with their activity predominantly restricted to Pseudomonas spp. Neither BD_H1 nor BD_H2 predators showed detectable lytic activity against K. pneumoniae or E. coli isolates tested.

To investigate potential genomic correlates of these phenotypic differences, we performed BLAST analyses of predation-associated genes (MAT adhesins, MIDAS, etc.) against assemblies from both BD_H1 and BD_H2 pools (Supplementary Table S7). All 12 core predation genes showed identical hits (100% identity/coverage) across both contigs, indicating no sequence variation in known predatory machinery. Core genome SNP analysis (Figure 2A) further confirmed minimal genomic divergence between isolates. These findings suggest that observed prey specificity differences likely arise from regulatory mechanisms, expression timing, or yet-unidentified genomic factors rather than variation in predation-associated coding sequences.

Variability between replicate experiments was observed for certain predator–prey combinations, likely reflecting the dynamic nature of predatory interactions and potential heterogeneity in prey surface structures. This variability underscores the importance of standardized assay conditions and replicate testing when evaluating predatory phenotypes (Caulton et al., 2024; Tyson et al., 2024).