Variation in genome content and predatory phenotypes between Bdellovibrio sp. NC01 isolated from soil and B. bacteriovorus type strain HD100

Defining phenotypic and associated genotypic variation among Bdellovibrio may further our understanding of how this genus attacks and kills different Gram-negative bacteria. We isolated Bdellovibrio sp. NC01 from soil. Analysis of 16S rRNA gene sequences and average amino acid identity showed that NC01 belongs to a different species than the type species bacteriovorus. By clustering amino acid sequences from completely sequenced Bdellovibrio and comparing the resulting orthologue groups to a previously published analysis, we defined a ‘core genome’ of 778 protein-coding genes and identified four protein-coding genes that appeared to be missing only in NC01. To determine how horizontal gene transfer (HGT) may have impacted NC01 genome evolution, we performed genome-wide comparisons of Bdellovibrio nucleotide sequences, which indicated that eight NC01 genomic regions were likely acquired by HGT. To investigate how genome variation may impact predation, we compared protein-coding gene content between NC01 and the B. bacteriovorus type strain HD100, focusing on genes implicated as important in successful killing of prey. Of these, NC01 is missing ten genes that may play roles in lytic activity during predation. Compared to HD100, NC01 kills fewer tested prey strains and kills Escherichia coli ML35 less efficiently. NC01 causes a smaller log reduction in ML35, after which the prey population recovers and the NC01 population decreases. In addition, NC01 forms turbid plaques on lawns of E. coli ML35, in contrast to clear plaques formed by HD100. Linking phenotypic variation in interactions between Bdellovibrio and Gram-negative bacteria with underlying Bdellovibrio genome variation is valuable for understanding the ecological significance of predatory bacteria and evaluating their effectiveness in clinical applications.

Introduction mechanisms governing predation in this strain (21). To investigate variation in genome content, 81 we compared the genome of NC01 to that of HD100 and other sequenced Bdellovibrio. 82 Currently, only a few completely sequenced predatory bacteria genomes are available in public 83 databases. In addition to HD100, there are only four other completely sequenced intraperiplasmic 84 Bdellovibrio genomes, including B. bacteriovorus W (22), which is the most phylogenetically 85 divergent from HD100 (23). The limited genome dataset highlights the importance of sequencing 86 additional Bdellovibrio, especially strains that are divergent from HD100. To define variation in 87 predatory phenotypes between NC01 and HD100, we determined prey range differences by 88 testing each strain against a panel of Gram-negative bacteria, and we measured predation 89 efficiency by quantifying viable E. coli ML35 during co-culture with each strain.

91
Isolation and identification of environmental bacteria for use as prey 92 We isolated bacteria from an urban freshwater stream in Providence, RI (lat 41.835, lon -  To identify these isolates, we sequenced the 16S rRNA gene following a previously 97 described procedure (24). Briefly, we used primers 63F (25) and 1378R (26), which amplify the 98 16S rRNA gene from many members of Domain Bacteria. Sanger sequencing of purified PCR 99 products was performed by GeneWiz (South Plainfield, NJ) using the same primers. We trimmed 100 and assembled reads using Phred/Phrap/Consed (27-29), then classified the assembled 16S 101 rRNA gene sequences using BLAST (30), the SILVA Incremental Aligner (31) and the in sterile HM buffer. Dilutions were then plated using a double agar overlay method. For each 127 dilution of the filtrate to be plated, we combined 100 µl of the dilution, 250 µl of Pseudomonas 128 sp. NC02 at 10 8 cfu/ml and 4.5 ml of molten HM top agar (0.6% agar). After vortexing, we 129 poured this mixture onto HM agar plates (1.5% agar). We allowed the top agar to solidify at 130 room temperature, then incubated plates at 28°C until we observed plaques in the lawn of 131 Pseudomonas sp. NC02. To begin the process of obtaining a pure isolate, we used a sterile 132 Pasteur pipette to extract single, well-defined plaques from the top agar. We combined each 133 plaque with 500 µl of sterile HM buffer, incubated the plaque suspension at room temperature 134 for at least 5 minutes, then vortexed it vigorously. 135 To culture predatory bacteria from the plaque suspensions, we made lysates combining 136 300 µl of each plaque suspension, 1.5 ml of an overnight culture of Serratia 0043 and 20 ml 137 sterile HM buffer. Serratia 0043 was also isolated from a soil environment, and switching to this 138 prey species during the isolation procedure selected for predatory isolates that were capable of 139 attacking both Pseudomonas and Serratia, which belong to different orders. We incubated 140 lysates at 28°C and 200 rpm. After 48 h, one lysate showed small, fast-moving cells, and we 141 used this lysate to begin two additional rounds of serial dilution, plating and plaque picking to 142 ensure a pure isolate. We used Serratia 0043 as prey during both of these rounds. After the final 143 round, we filtered the plaque suspension using a 0.45 µm filter, then combined 500 µl of filtrate 144 with 500 µl of sterile 50% glycerol and stored the stock at -80°C. 145 To identify NC01, we extracted the 16S rRNA gene sequence from the NC01 genome 146 based on the annotation generated by PGAP (see below for genome sequencing and annotation 147 methods). Using SINA 1.2.11 (31), we aligned this sequence against 16S rRNA gene sequences 148 extracted from nine genomes of obligate predatory bacteria belonging to Bacteriovorax, 149 Bdellovibrio Halobacteriovorax marinus BE01 CP017414, Halobacteriovorax marinus SJ NC_016620). We 154 then inferred a phylogenetic tree from the alignment with the online server RAxML Blackbox 155 (36), using 100 bootstrap replicates to estimate the confidence of the topology. 156 To further investigate the relationships among NC01 and other Bdellovibrio, we used an Negative stain electron microscopy of Bdellovibrio sp. NC01 161 To obtain samples of Bdellovibrio sp. NC01 for negative stain electron microscopy, we 162 added a small amount of the -80°C freezer stock to 20 ml of HM buffer mixed with 1.5 ml of an 163 overnight culture of E. coli ML35 grown in TSB. After 48 h of incubation at 28°C and 200 rpm, 164 we stained samples with 1% uranyl acetate (pH 4.2) following a previously described procedure 165 (24) and imaged the resulting specimens using a JEOL CX 2000 transmission electron 166 microscope. Ortholog groups were identified among all genomes based on the identification of bi-directional 211 best hits in the protein-coding genes using BLASTP for the pairwise sequence comparisons. The 212 pipeline requires an alignment of two amino acid sequences to have ³30% identity, ³70% 213 coverage, and an e-value of <0.001. To define the "core genome" of the seven Bdellovibrio 214 analyzed here, we identified the set of ortholog groups that each included at least one protein-215 coding gene from each genome. 216 We compared the core genome defined by our analysis with the set of 795 conserved 217 Bdellovibrio genes identified in (23). To do this, we extracted NCBI protein accession numbers associated with B. exovorus JSS from Additional File 4a, because JSS was the only genome for 219 which all 795 genes were listed. We then used R and bash scripts to determine whether these 220 accession numbers appeared in the ortholog groups identified as the core genome in our analysis.

221
To define the set of HD100 protein-coding genes implicated in predation by previous 222 studies, we extracted the HD100 locus tags listed in Supporting Information Figure S1 of (50) 223 and all HD100 locus tags with a value <0.5 in the column labelled Mini Tn-seq ECPL avg W in 224 Data Set S4 of (51). For each of these HD100 locus tags, we used R and bash scripts to 225 determine the corresponding GenBank protein ID (as listed in the NC_005363 annotation on 226 GenBank) and identify the ortholog group in our analysis that included this protein ID. We then 227 determined whether the ortholog group also included a NC01 protein-coding gene. If not, we 228 considered the gene missing in NC01.

229
For genome-wide pairwise comparisons of nucleotide sequences, we used the BLAST   To test prey range differences and plaque appearance, we plated 10-fold serial dilutions 241 of filtered lysate using a double agar overlay method as described above for the isolation 242 procedure. Each tested prey strain was grown in 120 ml TSB overnight at 28°C and 200 rpm, 243 with the exception of Escherichia strains, which were grown at 37°C. After centrifugation and 244 washing, we resuspended prey cell pellets in 5 ml of sterile HM buffer. Based on cfu/ml counts, 245 this yielded concentrations of at least 10 8 cfu/ml, which ensured a dense prey lawn for 246 visualizing plaques. For some prey strains, we added 500 µl of prey resuspension to top agar to 247 improve our ability to detect and observe plaques in the prey lawn. We incubated plates at 28°C 248 for 8 days before scoring plaque formation and performed three replicates for each prey strain. In  To test predation efficiency, we adapted a protocol described in (19). We centrifuged 252 filtered lysate at 8635 x g (8500 rpm in a Thermo Scientific A27 rotor) for 10 minutes. The  Figure 1). The average initial concentration estimated by direct cell counts was 1.05 x 10 7 279 cells/ml for NC01 and 2.3 x 10 7 cells/ml for HD100, whereas the average initial concentration 280 estimated by plaque counts (reported as plaque-forming units (pfu)/ml) at time zero was 1.17 x 281 10 7 pfu/ml for NC01 and 2.71 x 10 7 pfu/ml for HD100. These data show that direct cell counts 282 and plaque counts provide comparable results for quantification of Bdellovibrio.

284
Bioswale soil harbors predatory bacteria belonging to Bdellovibrio 285 Using Gram-negative soil bacteria isolated from a residential area as prey, we performed 286 an enrichment followed by three rounds of plaque picking to obtain a pure isolate of predatory 287 bacteria from a sample of soil collected at a bioswale on the Providence College campus. The 288 bioswale is a manmade landscape feature designed to collect and filter stormwater runoff from 289 nearby buildings. We designated the bioswale soil isolate as NC01. Using 1000x phase-contrast 290 microscopy of lysates combining NC01 and Gram-negative prey, we observed small, fast- and Bdellovibrio bacteriovorus HD100, which is the type strain of the species, showed 97% 302 similarity across 1,516 nucleotides. This value is above the 95% similarity cutoff typically used 303 to define bacterial genera (54), which supports identifying NC01 as Bdellovibrio. However, it is 304 below the 98.7% similarity cutoff typically used to define species (55), which suggests that 305 NC01 belongs to a different species than bacteriovorus. 306 To explore this further, we compared the average amino acid identity (AAI) of identified as B. bacteriovorus, cluster in a well-supported clade in the 16S rRNA phylogenetic 309 tree, and their amino acid sequences are highly similar (average AAI 99.2%). However, NC01 310 amino acid sequences are <70% similar on average to HD100 and 109J amino acid sequences, 311 supporting the conclusion that NC01 belongs to a different species than these two strains. NC01 312 amino acid sequences are even less similar on average to those of strain W, suggesting that 313 NC01 does not belong to the same species as this strain. Based on these analyses, we chose not 314 to designate a species for NC01, and we refer to this isolate as Bdellovibrio sp. NC01. The genome of Bdellovibrio sp. NC01 is larger than that of both HD100 and W (Table 2), 318 with a lower GC content than that of HD100, but slightly higher than that of W. We compared   SKB1291214, which is not complete (23). The core genome shared by these seven Bdellovibrio 331 strains included 795 protein-coding genes. We compared these to the PGAP core genome 332 described above and identified 778 protein-coding genes that were conserved in all eight 333 Bdellovibrio genomes. Of the 17 protein-coding genes that were identified as part of the core 334 genome by Oyedara and colleagues but not by our analysis, four were missing only in NC01.

335
The remaining 13 were missing in some of the other Bdellovibrio genomes according to our 336 analysis, likely due to differences in the parameters, such as alignment length, used by the two 337 software programs to cluster genes into ortholog groups.

338
To further investigate the four protein-coding genes identified as missing only in NC01, 339 we used BLAST to align representative amino acid sequences from the corresponding ortholog  In particular, BLASTP alignment showed 78% sequence similarity between the protein-348 coding gene annotated as a hypothetical protein (Bd1976/WP_011164428) and a 493 aa 349 sequence in NC01 that was annotated by both PGAP (DOE51_09405) and RAST (peg.1859).

350
The NC01 amino acid sequence is 161 aa longer than the representative amino acid sequence and 351 other sequences in the ortholog group due to a stop site difference. It was clustered in an ortholog 352 group by itself in the analyses of PGAP and RAST protein-coding genes, likely due to the 353 alignment length cutoff. Because the NC01 amino acid sequence shares high sequence similarity with members of the hypothetical protein ortholog group over two-thirds of its length, it may 355 have derived from the same ancestral gene. 356 Finally, the flagellar biosynthesis protein FliL was identified as missing only in NC01 357 because it was annotated as a pseudogene by PGAP (DOE51_16425). However, RAST 358 annotated it as a functional protein-coding gene (peg.3231). BLASTP alignment showed 89% 359 sequence similarity between a representative FliL sequence (Bd3329/WP_011165675) and the 360 NC01 amino acid sequence annotated by RAST, but over only 69% of the representative 361 sequence. This is due to differences in predicted start sites and therefore length between the 362 NC01 amino acid sequence (177 aa) and the representative amino acid sequence (236 aa). regions, TBLASTN did not produce any alignments to HD100 or W with >70% query coverage.

391
For those protein-coding sequences that did align to HD100 and/or W with >70% query 392 coverage, identity was low (<40%) for most of the sequences. Less than five protein-coding 393 sequences in each region aligned to HD100 or W at >70% coverage and >40% identity. This 394 suggests that the protein-coding sequences in the eight NC01 genomic regions are missing from 395 HD100 and W, instead of divergent.

396
In addition to comparisons against HD100 and W, we also used BLASTP to align amino 397 acid sequences annotated in each region by RAST against the non-redundant GenBank database.  Differential gene content between Bdellovibrio sp. NC01 and HD100 includes protein-428 coding genes that may be important for predation 429 In a pairwise comparison of protein-coding gene content between NC01 and HD100, we 430 found that ~32% of NC01 protein-coding genes are not found in HD100 and ~27% of HD100 431 protein-coding genes are not found in NC01 (Table 3). To assess broad functional classification 432 of these differentially present genes, we used RAST annotations of the genomes, primarily 433 because, in our experience, RAST predicts functions for more protein-coding genes, whereas 434 PGAP classifies more protein-coding genes as "hypothetical proteins". Additionally, both To assess variation in predatory phenotypes between Bdellovibrio sp. NC01 and B. 487 bacteriovorus HD100, we assayed prey range differences and predation efficiency. To determine 488 differences in prey range, we tested NC01 and HD100 against a panel of eight Gram-negative 489 bacteria isolated from different environments, including a freshwater stream and soil from a 490 residential area (Table 4). Plaque formation on a bacterial lawn indicates that the predatory 491 isolate is able to attack and lyse that particular Gram-negative prey strain. NC01 formed plaques on lawns of five of the eight prey strains tested, whereas HD100 formed plaques on lawns of all 493 eight prey strains. NC01 was not able to attack and lyse two prey strains belonging to 494 Acinetobacter and Raoultella that were isolated from a freshwater stream. In addition, although 495 NC01 formed plaques on E. coli ML35, it did not attack and lyse E. coli 0057.

496
To measure predation efficiency, we quantified viable E. coli ML35 by calculating cfu/ml 497 at different time points over 72 hours of co-culture with NC01 or HD100 (Figure 4) observed in all assay replicates, and cfu/ml plates showed no evidence of contaminant bacteria. 508 We quantified NC01 and HD100 after 72 hours of co-culture with E. coli ML35 to 509 determine if recovery of the ML35 population is associated with changes in predatory bacteria 510 population density. Using direct cell counts to estimate Bdellovibrio cells/ml at the start of co-511 culture and plaque counts to estimate pfu/ml after 72 hours of co-culture, we observed an 512 increase in the HD100 population of 1 log unit, whereas the NC01 population decreased by 1 log 513 unit ( Figure 5). When compared in validation tests, direct cell counts and plaque counts yielded 514 comparable results (Supplementary Figure 1), therefore the observed decrease in NC01 515 population is not due to a difference in methods for determining predator population density. 516 These data demonstrate that for both NC01 and HD100, active predatory bacteria are present 517 after 72 hours of co-culture with E. coli ML35, but, in contrast to HD100, the NC01 population 518 density decreases over the course of the predation efficiency assay.

519
To determine how predator population density during co-culture with E. coli ML35 520 compares to Bdellovibrio survival without prey, we quantified NC01 and HD100 over 72 hours  The only other species currently described within Bdellovibrio is exovorus; however, this species 547 uses an epibiotic predation strategy that differs from the cell invasion strategy observed by 548 microscopy for the soil Bdellovibrio isolate described here. Pending further characterization, we 549 designated the isolate as Bdellovibrio sp. NC01.

550
Because Bdellovibrio sp. NC01 is phylogenetically divergent from HD100 and other 551 intraperiplasmic Bdellovibrio with sequenced genomes, this isolate provides a valuable 552 opportunity to assess variation in genome content and predatory phenotypes among Bdellovibrio. 553 We clustered amino acid sequences from NC01, all five publicly available completely sequenced given Bdellovibrio genome are conserved among intraperiplasmic and epibiotic Bdellovibrio. 564 Analysis of the Bdellovibrio core genome also identified four protein-coding genes that 565 appeared to be missing in Bdellovibrio sp. NC01 but conserved in the other seven Bdellovibrio. 566 We confirmed the absence of two protein-coding genes using manual BLASTP alignments. One  Given these functional roles, the absence of the transporter permease in Bdellovibrio sp. NC01 574 may impact this strain's susceptibility to small inhibitory molecules.

575
For the remaining two protein-coding genes that appeared to be missing only in NC01, 576 manual BLAST alignments provided evidence of an ortholog encoded by the NC01 genome. In 577 particular, the PGAP annotation of NC01 flagged a gene encoding the flagellar biosynthesis 578 protein FliL as a pseudogene; therefore, it was not included in an ortholog group in our analysis. Bdellovibrio sequences in the database. This suggests that these Bdellovibrio strains may be 605 closely related to NC01, and the HGT events corresponding to acquisition of these six regions 606 occurred in the ancestral lineage before NC01 diverged from these strains. Addition of more 607 complete Bdellovibrio genomes to the database will improve our ability to reconstruct 608 evolutionary events such as horizontal gene transfer. in predatory phenotypes, it is essential to assign functions to protein-coding genes.

619
As the type strain, B. bacteriovorus HD100 has been the focus of multiple studies 620 investigating molecular mechanisms governing predation. We compared protein-coding gene 621 content between NC01 and HD100 and showed that ~30% of genes found in one of the strains 622 are missing in the other. These differentially present genes encode proteins from a wide range of 623 functional and structural families. To further investigate the possible impact of differentially 624 present genes on predation by NC01, we assembled a set of 343 HD100 protein-coding genes 625 identified as important for predation by previous studies using gene expression analysis or 626 mutagenesis (50, 51). About a quarter of these genes are missing in NC01. Most are annotated as 627 hypothetical proteins, but ten were annotated with functions that may be important for lytic 628 activity at different stages of the predatory lifecycle. In particular, four HD100 protein-coding 629 genes that are missing in NC01 are predicted to play a role in peptidoglycan metabolism.
Invasion of a prey cell by Bdellovibrio involves opening a hole in the prey cell peptidoglycan, 631 then resealing and reshaping peptidoglycan to form a rounded bdelloplast (11, 15). If NC01 is 632 missing enzymes that play a role in these processes, it may impact the ability of this strain to 633 invade prey cells and establish itself in the prey periplasm.

634
In addition to genome comparisons, we also compared predatory phenotypes between 635 Bdellovibrio sp. NC01 and B. bacteriovorus type strain HD100. Based on plaque formation on 636 lawns of eight different Gram-negative bacteria, NC01 has a more restricted prey range than 637 HD100. In particular, HD100 was able to attack and kill both tested strains of E. coli, whereas 638 NC01 was only able to attack and kill E. coli ML35. Other researchers have also observed that a 639 particular Bdellovibrio strain may vary in its ability to kill different members of the same species 640 (17, 63); however, the mechanisms governing these outcomes are not known. By demonstrating 641 that strain-level variation within a Gram-negative species determines susceptibility to predation,  reported as a characteristic of prey-independent mutants of Bdellovibrio; however, these plaques 663 are typically smaller than those formed by the wild-type strain and sometimes have a small 664 colony of host-independent Bdellovibrio at their center (64, 65). Neither of these features were 665 observed when comparing NC01 and HD100 plaques, suggesting that the NC01 plaque 666 phenotype is not the result of prey-independent growth.

667
The predatory phenotypes observed during NC01 co-culture with E. coli ML35 may be 668 the result of a mechanism specific to the prey, the predator, or a combination of the two. Other It is also possible that E. coli ML35 remains susceptible to predation throughout co-677 culture, but NC01 is deficient in some aspect of the predatory lifecycle. Our comparison of 678 NC01 and HD100 protein-coding gene content indicated that NC01 is missing some lytic 679 enzymes that may be important for predation. If these functions are not fulfilled by other protein- to prey, such as plastic phenotype resistance, and those specific to NC01 itself, such as deficient 687 or missing lytic enzymes, may not be mutually exclusive. Further work will test these 688 explanations.