Incomplete lytic cycle of a widespread Bacteroides bacteriophage leads to the formation of defective viral particles

Identification of B. thetaiotaomicron hankyphages with active DGR systems

To identify hankyphages integrated in Bacteroides genomes, we queried the genomic sequences of our previously described laboratory collection of 22 Bacteroides isolates (S1 Table [15],) with the reference hankyphage p00 from P. dorei HM179 [10]. Seven B. thetaiotaomicron strains out of the 22 (32%) contained hankyphage-like prophages, hereafter called HP+ strains, which were verified by PCR using hankyphage-specific primers (S1A Fig). Additionally, these strains were subjected to long read sequencing to obtain the complete prophage regions (S3 Table). One strain (jmh63) was discarded due to the presence of a fragmented hankyphage in the assembly (S1B Fig), resulting in six clinical lysogenic strains selected for further study. Their encoded hankyphages were generally conserved and mainly differed in the MuA-like transposase, putative tail fiber and DGR RT (Fig 1A). Additionally, all hankyphages grouped at the genus level (intergenomic similarity, IGS > 70%) with the reference hankyphage p00 (Fig 1A and S3 Table).

Fig 1. Characterization of Bacteroides thetaiotaomicron hankyphage laboratory collection.

(a) Genome alignment of our laboratory B. thetaiotaomicron hankyphages as well as the newly annotated reference Phocaelicola dorei HM719 hankyphage p00 using Clinker and Clustermap [16]. Gene annotations for genes that had a predicted function are indicated in black whilst hypothetical genes are numbered and annotated in gray. Ortholog coding sequences (CDS) are linked through bars colored by the percentage of amino acid identity. The DGR system components are highlighted in different colors (VR: Variable repeat (yellow), TR: Template repeat (blue), reverse transcriptase (pink), and Avd: accessory variability determinant (green)). (b) Boxplot depicting the nucleotide diversity (π) of the hankyphage VR region for individual overnight cultures of HP+ B. thetaiotaomicron strains with >30× coverage (reads spanning the complete VR region). We considered a significant DGR activity (red asterisk) if the mean diversity of positions targeted by the DGR (A in the TR) was greater than the mean of non-targeted positions +  2 sd [17]. (c) Bar plots indicating the proportion of A, C, T, and G nucleotides at each VR position, colored if they differ from the reference VR of each strain. Positions targeted by the DGR system are highlighted in dark blue. (d) Bar plot indicating absolute virion counts in the lysate of non-induced HP+ strains obtained by ddPCR. The HM719 strain is indicated in bold as it corresponds to the strain with the reference hankyphage p00. Results from three independent biological replicates, points indicate individual values and bars indicate their average, error bars indicate standard deviance from the mean, calculated using Prism. (e) Bar plot indicating absolute virion counts in B. thetaiotaomicron jmh43 lysates after the addition of potential inducers obtained by ddPCR. Results from 3 to 4 independent biological replicates, points indicate individual values and bars indicate their average, error bars indicate standard deviance from the mean, calculated using Prism. BHIS, indicated in bold, corresponds to the basal induction levels in culture medium. After confirming their normality, statistics were performed using a one-way ANOVA (mixed effect analysis with Geisser-Greenhouse correlation) and individual comparisons to the BHIS control were performed using Dunnet’s multiple comparison test. In the absence of indication, all comparisons showed non-significant differences. The individual quantitative values underlying this figure can be found in the S1 Data file also available at: https://figshare.com/s/3ff18cc2f6cc1edab0ae.

https://doi.org/10.1371/journal.pbio.3002787.g001

All hankyphage-like prophages encoded a predicted intact DGR system with the identification of identical TR sequences and diverse VRs (S3 Table) suggestive of DGR activity. To investigate whether hankyphage lysogenic strains grown in laboratory conditions underwent DGR mutagenesis, we analyzed sequencing data obtained from total DNA extracted from single overnight cultures. We then determined the nucleotide diversity (π) [18] in the VR region for the bacterial strains with at least 30 reads spanning the whole VR region (jmh42, jmh47, jmh51). We found that nucleotide diversity at positions targeted by the DGR, which corresponds to adenine residues in the TR region, was significantly higher than at non-targeted positions, with mutagenesis levels varying between strains (Fig 1B and 1C). This showed that the DGR diversifies the hankyphage tail fiber gene in our growth conditions.

Hankyphage virions are spontaneously produced by B. thetaiotaomicron isolates

Evidence of hankyphage DGR activity among our B. thetaiotaomicron strains prompted us to test their potential virion production. For this, the supernatants of stationary cultures were filtered, and viral DNA was purified. We then used multiplex digital droplet PCR (ddPCR) to obtain absolute quantification of the number of encapsidated hankyphage DNA molecules using a probe and primers located on the DGR RT. The absence of remaining genomic DNA was further checked by probing for the bacterial housekeeping gene rplB (S2 Fig and S2 Table). As a positive control, we used the lysate of P. dorei HM719 (carrying the p00 hankyphage), and as negative controls, we used the lysates of a B. thetaiotaomicron strain lacking the hankyphage (HP-) (jmh44) as well as the lysate of a genetically amenable B. thetaiotaomicron strain (jmh43) in which we deleted the hankyphage (jmh43∆HP). Using this method as a proxy for virion counting, we detected the presence of 10⁶–10⁸ virions/ml in the filtered supernatant of all HP+ B. thetaiotaomicron strains (Fig 1D), indicating that hankyphage virions were produced spontaneously in overnight cultures.

Previous studies reported a lack of plaque formation by hankyphage p00 on P. dorei naïve strains [10]. Consistently, none of the lysates of our collection of HP+ B. thetaiotaomicron strains yielded plaques on any of our HP-free strains nor on the jmh43 strain cured from its hankyphage (jmh43ΔHP). To test if we could increase virion production and maximize the possibility of new infection events, we focused on the genetically amenable B. thetaiotaomicron jmh43 strain, which carries a hankyphage highly similar to the hankyphage p00 (IGS =  99%) (Fig 1A). We tested a panel of potential inducers, including presence of salts (BHISP), different concentrations of the SOS-response inducer mitomycin C, various antibiotics (carbadox, erythromycin, ciprofloxacin), oxidative stress (H₂O₂ or oxygen exposure), different media imitating gut conditions (Gut Microbiota Medium (GMM) [19] or MiPro [20]), relevant gut cues (bile salts, lactose), osmotic shock (NaCl), iron chelators (EDTA, 2,2′-dipyridyl (BIP)), different pH (pH5 or pH9) and different temperatures (30 or 42 °C). None of these conditions revealed a significant increase in hankyphage virion production (Fig 1E), and the virions obtained did not lead to any plaque formation on naïve B. thetaiotaomicron strains (jmh43 ΔHP and jmh44). Overall, our results indicated that hankyphage virions were spontaneously produced by all tested HP+ strains from our collection and that none of the tested potential inducers increased virion production in B. thetaiotaomicron jmh43.

The CI-like repressor RepC_HP regulates hankyphage transcription

To increase virion production and raise the chances of observing infection events, we investigated if any genetic factors might trigger the hankyphage transition into the lytic cycle. For this, we looked for phage-encoded transcriptional regulators which could control the phage lysogeny. We noticed that one extremity of the hankyphage genome contained a small open reading frame, the location of which matched the position of the transposable Mu-bacteriophage repressor c (Repc) [21,22] (Fig 1B). This gene, therefore referred to as repc_HP, encoded a protein with a helix-turn-helix motif and shared a 37% amino acid identity to the λ phage CI repressor over the length of the full protein (Blastx). The Alphafold2 prediction of RepC_HP as a dimer showed a strong structural alignment with the crystal structure of the N-terminal DNA-binding domain of the λ phage CI dimer (root mean square of atomic positions, RMSD value =  2.44), suggesting that RepC_HP and the λ CI proteins could have similar DNA-binding properties (S3D Fig).

To test the involvement of repc_HP in the hankyphage virion regulation, we used the dead Cas9 (dCas9) technology [23–25] to inhibit repc_HP transcription in the HP + strain jmh43. We introduced a construct encoding an IPTG inducible dCas9 gene and a guide targeting the repc_hp coding region in a permissive site located outside of the hankyphage prophage within the jmh43 chromosome. In this strain (jmh43::dcas9-repc_HP), repc_HP is expected to be silenced upon IPTG addition. As a control, we cloned a non-targeting (nt) control containing a guide that does not target any B. thetaiotaomicron genes [23] at the same location (jmh43::dcas9-nt). To test the potential repressor function of repc_HP, we extracted and sequenced total RNA from exponential growing cultures of jmh43::dcas9-repc_HP and jmh43::dcas9-nt before IPTG induction (t =  0 h) and then at 1, 2, 6 and 22 h post-induction (hpi). Consistently with the constitutive virion production that we detected in the WT B. thetaiotaomicron jmh43 strain, we observed hankyphage transcripts in both strains already at t =  0 hpi (S3A Fig). Differential RNAseq analysis between the two strains at each time point showed that upon IPTG induction, dCas9 effectively silenced the expression of repc_HP, and led to a significant increase in the expression of most hankyphage genes from 1 h post-induction onwards (Fig 2A). Altogether, these results suggest that the product of repc_HP is a repressor of the hankyphage gene expression that maintains lysogeny.

Fig 2. The CI-like repressor RepC_HP regulates hankyphage transcription, replication and virion production.

(a) Log2 fold change in gene expression along the hankyphage genome in the jmh43::dcas9-repc_HP strain versus the jmh43::dcas9-nt control at different timepoints after IPTG-induced repc_HP silencing. Stars indicate significant differential expression, either upregulation (black star) or down-regulation (red star) calculated by Deseq2 and with a cutoff of >1 log2 fold change and a p-adjusted value b) Mapping of transposition events in jmh43::dcas9-repc_HP (red) and jmh43::dcas9-nt (blue) at 6 hpi (one biological replicate). (c) Quantification of transposition events of in jmh43::dcas9-repc_HP and jmh43::dcas9-nt at different timepoints (three biological replicates). (d) Quantification of virion production by ddPCR, CFU/mL (e) and absorbance (f) of jmh43::dcas9-repc_HP and jmh43::dcas9-nt* across different time points (three biological replicates). The individual quantitative values underlying this figure can be found in the S2 Data file also available at: https://figshare.com/s/3ff18cc2f6cc1edab0ae.

https://doi.org/10.1371/journal.pbio.3002787.g002

The RepC_HP repressor regulates hankyphage replication and virion production

Phage master repressors typically trigger the expression of phage genes, initiating a cascade of events contributing to the phage lytic cycle, including replication of the phage genome, virion assembly and exit from the host [22,26]. The hankyphage encodes a predicted MuA-like transposase and MuB-like ATPase [10] and it was hypothesized to replicate by transposition and insert randomly within the bacterial host genome, similarly to the Mu phage [10,27,28]. To test this, we investigated hankyphage replication following repc_HP silencing. Concretely, from the same jmh43::dcas9-repc_HP or jmh43::dcas9-nt cultures from which we extracted RNA, we also extracted whole-cell genomic DNA and examined the flanking chromosomal regions present at each end of the hankyphage genome. This revealed the presence of insertion events throughout the genomes of both strains and at a much higher frequency in the jmh43::dcas9-repc_HP strain compared to the jmh43::dcas9-nt control strain (Fig 2B), at every time point. These insertions culminated at an average of 2,504 transposition events at 22 hpi, compared to an average of 40 events in the control (Fig 2C). This formally demonstrated that the hankyphage could replicate by replicative transposition and that RepC_HP repressed phage replication.

To investigate if repc_HP silencing also had an impact on virion production, we performed a similar experiment whereby we extracted and quantified hankyphage virions by ddPCR from jmh43::dcas9-repc_HP and jmh43::dcas9-nt* (corrected strain, see material and methods section “Correction of jmh43::dcas9-nt strain”) growing cultures after IPTG induction. This showed an overall increase in the number of virions detected in jmh43::dcas9-repc_HP lysates compared to jmh43::dcas9-nt* , reaching an almost 1,000-fold increase at 22 hpi (Fig 2D). Silencing of repc_hp was accompanied by a significant decrease in colony-forming unit (CFU) counts compared to the control (Fig 2E), but a lack of effect on optical density (Fig 2F). This suggests that induction of the hankyphage lytic cycle represents a strong burden for the cells without leading to substantial cell lysis. To determine whether DNA packaged into hankyphage virions showed signs of DGR diversification, we extracted and sequenced DNA from DNase-treated PEG precipitated capsids of jmh43 and jmh43::dcas9-repc_HP lysates. We found DGR activity in both samples, with 7.93-fold (for jmh43::dcas9-repc_HP) and 3.32-fold (jmh43) more mutations at positions diversified by the DGR than non-targeted positions in the VR sequence (S4 Fig).

Overall, these results demonstrated that releasing RepC_HP repression leads to the induction of the hankyphage lytic cycle, including replicative transposition and virion production with active DGR mutagenesis.

Abnormal hankyphage virion assembly in jmh43 is partly recovered by repc
_HP silencing

To further examine the effects of repc_HP-silencing on hankyphage virions, we performed negative staining and transmission electron microscopy (TEM) on the lysates of jmh43 and jmh43::dcas9-repc_HP. TEM data presented in the initial description of the hankyphage p00 were suggestive of a tailed phage, consistent with the presence of predicted genes for neck and baseplate hub proteins [10]. Strikingly, we observed a total absence of tailed hankyphages in the lysate of jmh43 (Fig 3A). This lack of tails was also observed in the lysates of the reference P. dorei HM719 strain (S6 Fig). Only in the jmh43-dcas9-repc_HP lysate, three complete hankyphage particles could be observed from 10 different fields of view, as well as some loose tails. Both lysates also contained amorphous particles of different sizes ranging from 30 to 60 nm (S6 Fig) lacking the typical isometric appearance of mature capsids, which could be immature procapsids or vesicles (Fig 3A). Moreover, tailed hankyphages exhibited a greater capsid size than non-tailed particles (S7 Fig).

Fig 3. Abnormal hankyphage virion morphology and incomplete DNA packaging.

(a) Negatively stained TEM micrographs of lysates of jmh43 and jmh43::dcas9-repc_HP. Loose phage tails are indicated with black arrows. (b) Analytical ultracentrifugation sedimentation distribution of purified jmh43ΔHP and jmh43::dcas9-repc_HP lysates. The 280 nm spectrum has been normalized by the sum of the values per sample and is plotted on the left y-axis. The 260/280 nm ratio, indicative of the DNA/protein ratio is plotted on the right y-axis and values are indicated circles color-coded by sample as indicated by the legend. The figure depicts one out of two representative biological replicates. (c) Estimated sedimentation coefficients for the different expected assemblies in the jmh43 and jmh43::dcas9-repc_HP lysates. (d) Coverage map of DNA sequencing of PEG-precipitated and DNase-treated jmh43 and jmh43::dcas9-repc_HP lysates. Reads mapped to the hankyphage genome and normalized by the sum of mapped reads to that region. (e) Coverage map of DNA sequencing of genomic DNA samples of the jmh43::dcas9-repc_HP and jmh43::dcas9-nt strains at 22 hpi. Reads mapped to the hankyphage genome and normalized by the sum of mapped reads to that region. (f) Bar plot showing the ratio of hankyphage genome right-end and left-end copies quantified by ddPCR on ultracentrifuged capsids of HP+ strains. The individual quantitative values underlying this figure can be found in the S3 Data file also available at: https://figshare.com/s/3ff18cc2f6cc1edab0ae.

https://doi.org/10.1371/journal.pbio.3002787.g003

To confirm the nature of the different macromolecules present in the repc_HP-silenced lysate, we performed analytical ultracentrifugation (AUC) and analyzed the sedimentation profiles of the concentrated lysates of jmh43::dcas9-repc_HP and the negative control jmh43ΔHP (Fig 3B). To decipher the type of assembly in each peak, the expected sedimentation values of different hankyphage assemblies were modeled using the dimensions observed in the negative stain micrographs (Fig 3C) (See methods section ‘Analytical ultracentrifugation’). For each peak, the 260/280 ratio was calculated to detect the presence or absence of DNA (Fig 3B). Different assemblies could be detected in the jmh43::dcas9-repc_HP sample. One sedimented at ~ 62S with a high protein content (low 260/280 ratio), consistent with the expected sedimentation value of hankyphage tails. Another sedimented at 120S with a high protein content, consistent with empty capsids. Two broad peaks could be seen from 210S to 320S and then from 400S to 520S, both with a high DNA content, and which could correspond to procapsids or tail-less capsids, respectively (Fig 3B and 3C). The low abundance of fully assembled virions seen in our microscopy data makes it unlikely that we can detect them with AUC. These results agreed with the observation of detached tails and putative immature capsids in the jmh43::dcas9-repc_HP negative staining images. The control jmh43ΔHP sample did not show any high DNA content particles; instead, it showed different assemblies with a low (dcas9-repc_HP.

The jmh43 strain encodes four additional predicted prophages, which can be potentially induced and confound results (S5 Table). Our RNAseq data detected expression for only some genes of these prophages and not of the full genome – in contrast to the hankyphage where all genes were expressed (S7 Fig). In line with this, mass spectrometry analysis of the jmh43::dcas9-repC_HP lysate revealed the presence of most of the hankyphage proteins and an absence of proteins from the other four predicted prophages (S6 Table), suggesting that only hankyphage virions contribute to our observations. Additionally, we did not find an enrichment of coverage of other chromosomal regions beyond the hankyphage loci in the DNA extracted from the capsids for both jmh43 and jmh43::dcas9-repc_HP strains.

Given the presence of some tailed virions in jmh43::dcas9-repc_HP lysates, we decided to test their infectivity. However, plaque assays using jmh43::dcas9-repc_HP lysates yielded no plaques on our HP-free strains nor jmh43ΔHP. Considering the role played by the B. thetaiotaomicron capsule for phage infection, we hypothesized that specific capsule types that might only be transiently expressed might allow hankyphage adsorption [5]. B. thetaiotaomicron contains 8 different capsule loci that are regulated in an on/off manner and generate extensive surface diversity [29]. We tested the infection potential of jmh43::dcas9-repc_HP lysates as well as jmh43ΔHP lysates as a negative control, on a collection of capsule mutants of the hankyphage-free B. thetaiotaomicron VPI5482 strain, either expressing one of the eight single capsule types (cps1–8), none at all (△cps1–8), or heterogeneously expressing the eight (VPI548) [30]. We did not find any significant reduction in optical density for the non-capsulated mutant exposed to the jmh43::dcas9-repc_HP lysate compared to the jmh43ΔHP lysate control (S8A Fig). Neither did we observe plaque formation on any of the mutants (S8B Fig).

Given the temperate nature of hankyphages, it could be possible that frequent and rapid lysogenization events impaired subsequent infections. We thus tested whether jmh43 and jmh43::dcas9-repc_HP virions could lysogenize new strains. For this, we inserted an erythromycin-resistance encoding gene (ermG) with its own promoter on the repc_HP end of the jmh43 hankyphage, downstream of repc_HP in the WT strain (jmh43-HP-Erm) as well as the repc_HP silenced strain (jmh43::dcas9-repc_HP-HP-Erm) (S8C Fig). We incubated the filtered supernatant of an overnight culture of strain jmh43-HP-ermG, containing ermG-tagged hankyphage, with HP-free strains (jmh∆43∆HP and jmh44) and the VPI5482 capsule mutants for 2 h followed by plating in the presence of erythromycin. PCR verification revealed that the few colonies obtained did not contain the ermG gene (S8D Fig), suggesting they probably were spontaneous HP-free erythromycin resistant mutants and indicating an absence of any new hankyphage lysogenization events. Additionally, to better mimic the gut’s physiological conditions which might facilitate phage infection, we performed the same transduction experiment using mouse fecal extracts, but this did not yield any transductants either.

Overall, our results show that the B. thetaiotaomicron jmh43 hankyphage spontaneously produces defective virions. However, despite the presence of complete virions after repc_HP-silencing, jmh43 hankyphages are non-infective in our laboratory conditions.

Defects in hankyphage genome packaging in vitro

DNA sequencing of jmh43 and jmh43::dcas9-repc_HP capsids revealed an enriched coverage in the entire hankyphage region, from which we could determine that this transposable phage packages 29–31 bp of the host chromosome on its 5′ end, and a variable fragment of 1–2.3 kb on its 3′ end, consistently with a headful packaging mechanism (S9 Fig).

Surprisingly, we observed that the hankyphage genome was abnormally encapsidated in the two samples, with the left, repc_HP-end (5′) of the hankyphage genome covered 40-fold times more than the right, methyltransferase (3′) end, with the coverage dropping gradually (Fig 3D). This contrasted with the uniform coverage obtained when sequencing total DNA extracted at different time points from the jmh43::dcas9-repc_HP or jmh43 cultures (Fig 3E). Overall, our data shows that the jmh43 hankyphage DNA packaging is initiated at the left end of the phage but is randomly interrupted in the first 15 kb or ejected from the capsid, with only a small fraction of capsids being fully packaged. To investigate whether this phenomenon was specific to the hankyphage of strain jmh43 or a more general feature of hankyphages, we quantified the abundance of DNA at the left (MuA-like gene) and right end (DGR RT gene) of the hankyphage in virions present in the PEG-precipitated capsids of other HP+ B. thetaiotaomicron strains. In all the strains tested, the left end/right end ratio was even greater than the one detected in jmh43 (S10 Fig). To check that these differences were not due to virion damage during PEG concentration, we performed the same analysis after phage ultracentrifugation and found the same trend, albeit to a lesser extent (Fig 3F). Altogether, these results suggest that the hankyphages in all HP+ B. thetaiotaomicron strains of our collection have a packaging defect or are unstable in laboratory conditions.

Genomic analysis of hankyphages suggests active mobility

Our results could suggest that at least some B. thetaiotaomicron and P. dorei hankyphages are no longer able to transfer horizontally due to structural defects. These defects could be the consequence of the burden potentially associated to the hankyphage spontaneous induction and transposition. However, deletion of the hankyphage generally did not significantly increase the growth of the host compared to the WT strains (S11 Fig), showing that the hankyphage does not pose a strong burden to the bacteria.

To gain a comprehensive understanding of the hankyphage distribution and potential mobility within Bacteroides and Phocaeicola, we compiled a dataset of 136 complete genomes from the NCBI, supplemented with the 25 strains from our laboratory collection (S4 Data). Out of these 161 strains, we found 26 (16%) encoding hankyphage-like prophages, with 20 of them grouping at the genus level (IGS > 70%) and 9 to the species level (IGS > 95%) with the reference hankyphage p00 [10]. Hankyphages closely related to this reference genome (at least 75% nucleotide identity, with most hankyphages having >95% nucleotide identity) were scattered across the phylogeny of Bacteroides and Phocaeicola (Fig 4). We then compared the prophage and bacterial phylogenies to test for potential evolutionary history correlation, and we observed that the phage phylogeny was not concordant with the phylogeny of the hosts and had multiple crossovers (Fig 4). Together with the presence of very similar hankyphage genomes in evolutionary distant bacterial hosts (S12 Fig), this suggests that there were recent events of hankyphage horizontal transmission. Moreover, analysis of the coverage pattern of hankyphages in metaviromes [10] did not reveal a drop in coverage toward the 3′ end of the hankyphage genome (S13 Fig). This suggested that the defect in capsid DNA packaging that we observe in our laboratory conditions does not occur in vivo. Overall, these results support the hypothesis that hankyphages are active and can transfer horizontally.

Fig 4. Phylogenetic analysis of Bacteroides hankyphages.

Tanglegram showing Bacteroides (left) and hankyphage phylogenies (right) constructed from their respective core pangenomes [31], as well as crossover events representing the presence of Hankypages integrated in the genomes of the respective bacterial hosts. Hankyphage presence (identity above 75%, relative to the reference Phocaelicola dorei HM719 hankyphage p00) is indicated with a blue box, with darker blue boxes indicating a coverage of at least 60%, and lighter blue boxes indicating a coverage of less than 60% of the reference hankyphage genome. Tree files for Bacteroides and hankyphages and their respective IDs can be found in supplementary S17, S18 and S4 Data, respectively, available at: https://figshare.com/s/3ff18cc2f6cc1edab0ae.

https://doi.org/10.1371/journal.pbio.3002787.g004

Bacterial strains and growth conditions

Bacteroides strains were grown in brain-heart infusion medium supplemented (BHIS) with cysteine 0.1%, hemin 5 μg/mL and sodium bicarbonate 0.2% (NaCOH₃). Cultures were incubated statically at 37 °C in a C400M Ruskinn anaerobic chamber (95% N₂, 5% CO₂) or in anaerobic bags (GenBag anaero, Biomerieux). E. coli strains were grown in Luria-Bertani (LB) broth (Corning) supplemented with or without ampicilin 100 μg/mL and incubated at 37 °C aerobically with 180 rpm shaking.

Construction of
B. thetaiotaomicron
marker-less deletion mutants

Deletion mutants of B. thetaiotaomicron were performed by allelic exchange using the pLGB13 vector [52]. Briefly, 500 base pairs upstream and downstream regions of the target sequence were cloned into the pLGB13 vector by Gibson assembly and transformed into E. coli S17 λpir. This strain was then used to deliver the vector to B. thetaiotaomicron by conjugation. Matings were carried out by mixing exponentially grown cultures of the donor and the recipient strains at a 2:1 ratio and placing the mixture on BHIS agar plates at 37 °C under aerobic conditions overnight. The next morning, bacteria were plated on selective BHIS agar supplemented with erythromycin (15 µg/mL) for selection of B. thetaiotaomicron transconjugants that underwent the first recombination event and gentamicin (200 µg/mL) to ensure the exclusion of E. coli growth. Erythromycin resistant colonies of B. thetaiotaomicron were then grown overnight without antibiotics to promote plasmid loss and subsequently subjected to a second round of selection by plating on BHIS agar plates supplemented with anhydrotetracycline (for selection of double recombined colonies). The resulting deletion mutants were confirmed by PCR using external primers and by sanger or nanopore amplicon sequencing, or whole-genome illumina sequencing.

dCas9-mediated silencing

An easily modifiable vector based on the pMM704 plasmid [23] was constructed by Gibson assembly where all BsaI sites were removed and an mCherry gene flanked by bsaI sites was introduced at the position of the guide insertion. Oligonucleotides were then designed to allow cloning of the desired guide between two BsaI sites on the vector. Oligonucleotides were phosphorylated and annealed to reprogram the guide carried by the vector using the Golden Gate assembly method. The mixture was then dialized and introduced into E. coli S17 λpir by electrophoration. The clones were confirmed by PCR and nanopore sequencing and if correct, were conjugated into B. thetaiotaomicron as described above. Bacteria were then selected on selective BHIS agar supplemented with erythromycin (15 µg/mL) and gentamicin (200 µg/mL) to ensure the exclusion of E. coli growth and the resulting mutants were confirmed by PCR using external primers and sanger or nanopore amplicon sequencing, or whole-genome illumina sequencing.

Correction of jmh43::dcas9-nt strain

Differential expression analysis on the RNAseq data revealed a differential over-expression of the dcas9 gene in jmh43::dcas9-repc_HP compared to the jmh43::dcas9-nt strain (S3C Fig), which led us to find a 6-bp deletion in one of the lac operators of the IPTG-inducible promoter controlling the expression of dCas9 in jmh43::dcas9-nt. This strain was re-constructed to correct the promoter sequence (jmh43::dcas9-nt*) (S4 Table) and further experiments were performed with both versions of the strain, yielding similar results.

Hankyphage induction

Overnight cultures were diluted to an absorbance of OD₆₀₀ = 0.05 and grown anaerobically in 5 mL cell culture tubes at 37 °C until exponential phase. For time point experiments, cultures were grown in 20 mL volume in 50 mL falcon tubes. Inducers were added to the following concentrations: BHISP (standing for BHIS-Phage media, contained the following salts: CaCl₂ 1 mM and MgSO₄ 10 mM), mitomycin C 2 μg/mL or 5 μg/mL, carbadox 8 μg/mL, erythromycin 5 μg/mL, ciprofloxacin 0.1 μg/mL, bile salts 0.5%, NaCl 0.5 M, EDTA 0.2 mM, bip (2,2′-dipyridyl) 0.2 mM, hydrogen peroxide (H₂O₂) 0.5 mM, lactose 0.4%, 42, 30 °C. For the cultures induced in GMM or MiPro medium, overnight cultures were adjusted to an OD600 = 0.3 and then washed and resuspended in either medium supplemented with 0.4% glucose. The cultures were then incubated anaerobically at 37 °C for 18 h. The cultures were then centrifuged at 7142g for 15 min at 4 °C, the supernatant was then filtered using 0.45 μm filters and the lysate stored at 4 °C.

Lysis curves

Overnight cultures were diluted to an absorbance of OD₆₀₀ = 0.05 and grown anaerobically in 5 mL cell culture tubes at 37 °C for 2 h until exponential phase. Lysates of different strains were mixed in a 1:1 ratio in a total volume of 200 µl in 96-well Greiner flat-bottom plates with the exponentially-grown bacterium inside the anaerobic station, and then sealed with an adhesive sealing sheet (Thermo Scientific, AB0558). The plates were then incubated in a TECAN Infinite M200 Pro spectrophotometer for 24 h at 37 °C. Absorbance at OD₆₀₀ was measured every 15 or 30 min, after a 200-second orbital shaking with a 2 mm amplitude.

Transductions

Overnight cultures were diluted to an OD600 = 0.05 and grown anaerobically in 5 mL cell culture tubes at 37 °C for 2 h until exponential phase. One mL of cells were then pelleted and were either resuspended with 200 µl of the phage lysate or in a 1:1 ratio mixture of the phage lysate and mouse faecal filtrate. These were then incubated anaerobically and statically at 37 °C for 2 h. The mixture was then plated on erythromycin-supplemented BHIS plates and grown for 48 °C at 37 °C anaerobically. The colonies obtained were re-streaked on erythromycin-supplemented BHIS and single colonies were tested by PCR for the presence of the hankyphage and the ermG gene.

CFU counts

Cultures were diluted serially in 10-fold or 100-fold dilutions in BHIS and 5 µl drops were placed on BHIS or BHIS-erythromycin solid media and grown for 48 h at 37 °C anaerobically. Colonies were counted from the most serially diluted spot containing >3 and

Growth curves

In Greiner flat-bottom 96-well plates, overnight cultures were diluted to a final absorbance of 0.1 (OD₆₀₀) within 150 µ L BHIS in anaerobic conditions and a plastic adhesive film (adhesive sealing sheet; Thermo Scientific, AB0558) was used to seal the plate. The plates were then incubated in a TECAN Infinite M200 Pro spectrophotometer for 24 h at 37 °C. OD₆₀₀ was measured every 15 min, after a 200-s and 2 mm amplitude orbital shaking.

Virion quantification by digital PCR (ddPCR)

Phage lysates were DNase treated to remove remains of genomic DNA by adding 2 μL of TURBO DNase (Thermo Fisher, REF #AM2238) and 10 μL of 10× TURBO Dnase Buffer to 88 μL of lysate, incubating at 37 °C for 1 h and then heat-inactivating the enzyme and opening the phage capsids by a 10-min 99 °C treatment. ddPCRs were performed as triplex reactions. For quantification of the virion production of HP+ strains, as well as the quantification of hankyphage right-end and left-end copies, we used a probe and primers specifically targeting a conserved region of the hankyphage muA gene and another targeting the rt gene, as well as the rplB-targeting genomic control (S2 Table). Reactions on jmh43 lysates obtained from different inducers as well as on jmh43::dcas9-repc_HP and jmh43::dcas9-nt were performed with a probe and primers specific to the hankyphage rt gene, a second probe specific for the housekeeping gene rplB, and a third probe specific for the 3′ flanking region of hankyphage and the jmh43 chromosome (S2 Table).

Reactions were prepared either using the Naica PCR mix (R10056) or the Perfecta Multiplex qPCR ToughMix (Quanta Biosciences, Gaithersburg, MD, USA) and placed on either Sapphire or Ruby chips (Stilla technologies). Digital PCR was conducted on a Naica Geode using the following steps: droplet partition (40 °C, atmospheric pressure AP to +  950 mbar, 12 min), initial denaturation (95 °C, + 950 mbar, for 2 min), followed by 45 cycles at (95 °C for 10 s and 58 °C for 30 s), droplet release (down 25 °C, down to AP, 33 min). Image acquisition was performed using the Naica Prism6 reader. Images were then analyzed using the Crystal Reader (total droplet enumeration and droplet quality control) and the Crystal Miner software (extracted fluorescence values for each droplet) (Stilla Technologies). Values were then transformed from copies/µl to copies/mL.

Identification and annotation of hankyphages

All complete Bacteroides and Phocaeicola genomes from the NCBI were obtained with ncbi-genome-download (accession: 08/08/23) [53]. Additionally, we included 22 Bacteroides isolates available in our laboratory [15] and four P. dorei strains obtained from the Biodenfense and Emerging Infections (BEI) via the American Type Culture Collection (ATCC), resulting in a total of 161 genome sequences (S4 Data). The classify workflow from GTDB-Tk v2.1.1 [54] was used to correct the taxonomy classification of the bacterial genomes. geNomad [55] was used to detect prophages which were then queried against the hankyphage reference genome (BK010646, 42,831 bp) using blastn [56] (-evalue 1e-5 -max_hsps 1). Candidate hankyphage-like prophages were extracted from hits with a minimum alignment length of 10% with the reference phage with seqtk (https://github.com/lh3/seqtk). A blastn search independent of geNomad identification revealed some additional HP+ clinical strains. In these strains, the phage was highly fragmented in the assembly, which prevented its detection by geNomad in the initial analysis. HP + laboratory strains were subjected to nanopore sequencing to improve the quality of the genome assembly (ArrayExpress accession number: E-MTAB-14724). Hybrid assemblies of Illumina and nanopore reads were obtained with Unicycler v. 0.4.8 [57,58]. Hankyphage-like prophages were extracted as above and further annotated with pharokka [58] (-p prodigal) and phold, an annotation pipeline that uses protein structures to annotate phage genes [59,60]. Genes with low confidence were manually curated using Foldseek [61]. Annotated hankyphages can be found at the NCBI database, accession: PRJNA1200653. Extracted hankyphage-like prophages were further compared with VIRIDIC which calculates IGS to classify phages at the genus and species level [62]. Clinker was used to obtain orthologous genes and comparison between hankyphage-like genomes [16]. Defense systems were identified with DefenseFinder [44,63,64].

Alphafold2 modeling of RepC_HP

The complete amino acid sequence of RepC_HP from the jmh43 hankyphage was submitted to the Alphafold2 collab server as a dimer with the following parameters: num_relax = 0, template_mode=none. The model with highest confidence was aligned to the N-terminus fragment of a dimeric crystal structure of the λ CI repressor (PDB3kz3) [65] using the PyMOL function “align”.

RNA isolation and sequencing

Overnight cultures were mixed with RNAprotect (Qiagen) in the anaerobic chamber to prevent RNA degradation and bacteria were lysed using QIAGEN Proteinase K and TE buffer containing lysozyme. Total RNA was extracted using the Direct Zol kit (Zymo Cat. R2050) according to the manufacturer’s instructions and treated with Dnase I from the same kit. RNA concentration, quality, and integrity from three independent replicates were checked using the 4195 Tapestation system (Agilent).

RNAseq analysis

Sequencing was performed by the Biomics platform at the Institut Pasteur. Ribosomal RNA depletion was performed using the Bacteria RiboZero Plus kit (Illumina) with a Bacteroides thetaiotaomicron custom ribodepletion probe. From rRNA-depleted RNA, directional libraries were prepared using the Illumina Stranded total RNA Prep ligation preparation kit following the manufacturer’s instructions (Illumina). Libraries were checked for quality on Fragment Analyzer (Agilent). Quantification was performed with the fluorescent-based quantitation Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). 750pM of pooled library were sequenced in pair-end 100 cycles on a NextSeq2000 Illumina sequencer. The RNA-seq analysis was performed with Sequana 0.14.1 [66]. We used the RNA-seq pipeline 0.15.2 (https://github.com/sequana/sequana_rnaseq) built on top of Snakemake 7.8.5 [67]. Briefly, reads were trimmed from adapters using Fastp 0.20.1 [68] then mapped to the jmh43::dcas9-repc_HP genome (extracted at time point 0 hpi) using bowtie2 v.2.4.2 [69]. The reference jmh43::dcas9-repc_HP-t0hpi genome (accession: ERR14225767) was annotated using Prokka [70] and the hankyphage hypothetical genes were attributed specific names based on their coordinates. FeatureCounts 2.0.1 [71] was used to produce the count matrix, assigning reads to features using corresponding genome annotation with strand-specificity information. Quality control statistics were summarized using MultiQC 1.11 [72]. Statistical analysis on the count matrix was performed to identify differentially regulated genes. Clustering of transcriptomic profiles were assessed using a Principal Component Analysis. Differential expression testing was conducted using DESeq2 library 1.30.0 [73] scripts indicating the significance (Benjamini-Hochberg adjusted p-values, false discovery rate FDR p-adjusted value S2 Data and S6 Data respectively. Raw reads can be accessed at the EMBL functional genomics data collection (ArrayExpress) (accession: E-MTAB-14274).

Genomic DNA sequencing

Whole-cell genomic DNA was extracted from 2 mL of pelleted cultures using the Qiagen Blood and Tissue Kit (REF #69504) following the instructions for Gram negative bacteria. Samples were then RNAse-treated and resuspended in 50 µl of DNase and RNase-free water. The samples were then sequenced at the Institut Pasteur Biomics platform using a Nextseq2000 Illumina sequencer with a P1 flowcell. Libraries were prepared using the Illumina Truseq DNA PCR-Free kit. Sequencing was performed as paired-end in 300 cycles.

Quantification and mapping of transposition events

Raw reads from genomic DNA (available at ENA—accession: PRJEB85302) were mapped to the jmh43 hankyphage prophage using bwa mem v. 0.7.17-r1188 [74], following delimitation of the prophage ends using capsid data. If the prophage had integrated into a different chromosomal region, reads at the ends would be clipped. Soft-clipped reads were extracted with the Pysam software package [75]. The clipped regions (≥ 15 nt) were subsequently searched into the jmh43 bacterial chromosome using blastn (-task blastn-short -qcov_hsp_perc 100 -perc_identity 100). Clipped reads that were mapped into multiple genome locations were discarded (76].

Mutations at VR region

DGR components (RT, TR and VR regions) were detected as previously from genome assemblies [77]. Briefly, hmmsearch v. 3.4 was used to detect RT based on matches to HMM profiles of curated RT. Then, repeats around the candidate RT were detected with blastn (option -word_size 8 -dust no -gapopen 6 -gapextend 2). Finally, the selection of the TR and target genes were based on the repeat patterns (e.g., mutation bias and distance to the RT). The pipeline used is available at https://bitbucket.org/srouxjgi/dgr_scripts/src/master/DGR_identification/. The target region, which includes the VR, was aligned with mafft v.7 [78]. The resulting nucleotide alignment was used to design conserved primers flanking the VR region of the clinical isolates (F: “ATATCCACTGCCTCC”, R: “GCCAGACGAGAACCG”). Primers were used to specifically collect reads spanning the complete VR region using SeqKit v.2.3 amplicon function [79]. The resulting amplicons were aligned with mafft, and the genotypes analyzed with R. To distinguish between DGR activity and background mutations/errors, we computed nucleotide diversity at each position as π =  1 − (A² + T² + C² + G²). We considered a VR to be diversified in a given strain if mean(π,A) > mean(π,nonA) + 2 * sd(π,nonA) [17]. Alternatively, we used a mapping approach to include also reads that span some region of the VR. Briefly, reads matching VR flanking primers were first filtered using bbduk.sh [80] (bbduk options k = 15 mm = f). The resulting reads were mapped against the VR region ± 100 nt using bbmap to obtain a global alignment of the reads instead of the soft-clipped alignment of local alignments (bbmap options “vslow minid=0”) [81]. A mpileup file was obtained using samtools mpileup (options -f and -A to annotate reference alleles and keep anomalous read pairs, respectively). A custom R script was then used to obtain allelic frequencies and nucleotide diversity from the raw mpileup file.

Transmission Electron Microscopy (TEM)

All micrographs shown in the main text were prepared by mixing the lysates with paraformaldehyde (4% final concentration) for 24 hours and then incubating 4 µl of lysate for 1 min on CF300 carbon-coated grids (previously glow discharged 2 mA, 1 min in ELMO). The grids were rinsed four times in MilliQ water drops and contrasted using 2% uranyl acetate. Images were acquired with a Tecnai (ThermoFisher) T12 operating at 120 kV and a RIO 13 camera (Ametek).

For the micrographs in S5 Fig and S15 Fig, 4 µl of lysate was incubated for 1 min on CF300 carbon-coated grids (previously glow discharged 2 mA, 1 min in ELMO). The grid was rinsed four times in MilliQ water drops and contrasted using 2% uranyl acetate. Images were acquired with a Tecnai (ThermoFisher) T12 operating at 120 kv and a RIO 13 camera (Ametek). Particle size measurements were performed using FIJI [82].

Analytical ultracentrifugation (AUC)

Overnight cultures were diluted to an absorbance of OD600 = 0.05 and grown anaerobically in 15 mL cell culture tubes at 37 °C until exponential phase. The cultures were then incubated anaerobically at 37 °C for 18 h with or without addition of IPTG (500µM). The cultures were then centrifuged at 7,142g for 15 min at 4 °C, the supernatant was then filtered using 0.45 μm filters. Lysates were then concentrated, and medium exchanged for PBS (no calcium nor magnesium) by placing the lysate in amikon 100 kDa cutoff columns and spinning for 5 m at 3,000g at 4 °C or until the sample was concentrated to 500 µl twice, keeping the top fraction. 130 µl of concentrated lysate was mixed with DNase buffer (1× final concentration) and treated with DNase Turbo (3 µl) for 2 h at 37 °C. The DNase treated lysates were then resuspended in 1 mL of PBS and spun down in new 15 mL Amikon 100 KDa columns to a final 500 µl volume to remove the DNase enzyme. Samples were centrifuged at 6,000 rpm in an Optima AUC analytical ultracentrifuge (Beckman Coulter), at 20 °C in an eight-hole AN 50–Ti rotor equipped with 12-mm double-sector aluminum epoxy centerpieces. Detection of the biomolecule concentration as a function of radial position and time was performed by absorbance measurements at 260 nm, 280 nm and by interference detection. Ultracentrifugation experiments were performed in PBS. Sedimentation velocity data analysis was performed by continuous size distribution analysis c(s) with invariant D [83] using Sedfit 16.36 software. All the c(s) distributions were calculated with a fitted diffusion and a maximum entropy regularization procedure with a confidence level of 0.68. To evaluate the amount of DNA 260/280 absorbance ratio were calculated at the top of each peak; to note ratio below 1 shows the presence of less than 5% of DNA sample. Therefore, for ratio below 1 the amount of DNA is negligeable [84]. Buffer viscosity and density were calculated from the Sednterp software. Modeling of the different assemblies was evaluated with dimensions measured by TEM using spherical modeling as suggested by previous literature [85,86].

Mass spectrometry

PEG-precipitated capsids (see ‘capsid DNA extraction and sequencing’ section) from one biological replicate of a jmh43::dcas9-repc_HP lysate were mixed with 1× Laemli buffer (Bio-Rad #1610747) and boiled at 100 °C for 10 min. After separation by denaturing SDS–PAGE using precast TGX 4%–15% gradient gels (Bio-Rad), the gel bands of interest were cut and washed using acetonitrile (ACN)/water-based dehydration/rehydration steps. Disulfide reduction was performed by adding 200 μL of a solution containing 10 mM dithiothreitol in ammonium bicarbonate (AMBIC) 50 mM to dehydrated gel bands for 30 min at 56 °C. Thiol alkylation was performed by incubating the previous samples in 20 μL of a 500 mM iodoacetamide solution in AMBIC for 30 min at room temperature. Reagents were removed using the previous washing procedure. Proteins were digested with a solution of 0.05 μg/μL of trypsin/Lys-C (Promega) in 100 μL of AMBIC overnight at 37 °C. The resulting peptide mixture was filtered and acidified with trifluoroacetic acid at a final concentration of 0.1% (v/v). For each sample, 6 μL were loaded on a C18 cartridge (Dionex Acclaim PepMap100, 5 μm, 300 μm i.d. × 5 mm) and eluted on a capillary reversed-phase column (Waters nanoEaze M/Z Peptide CSH C18, 1.8 μm, 75 μm i.d. × 25 cm) at 220 nL/min, with a gradient of 2% to 38% of buffer B in 60 min (buffer A: 0.1% aq. Formic Acid/ACN 98:2 (v/v); buffer B: 0.1% aq. Formic Acid/ACN 10:90 (v/v)), using a nanoRSLC U3000 system coupled with a Tribrid Orbitrap Eclipse mass spectrometer (ThermoFisher Scientific) using a “Top Speed” data-dependent acquisition MS experiment: 1 survey MS scan (400−2,000 m/z; resolution 70,000) followed by the highest possible number of MS/MS scans within 3 s, in the linear ion trap (dynamic exclusion of 30 s). Protein identification was performed with Proteome Discoverer 2.4 using SEQUEST against the jmh43 genome, with the following parameters: methionine oxidation, asparagine/glutamine deamidation as variable modifications, cysteine carbamidomethylation as fixed modification, MS1 error tolerance set at 10 ppm and MS/MS error tolerance at 0.6 Da, protein identification at FDR 1%. The identified proteins were compared against the proteome of the jmh43 strain using blastp (-e-value 1E−4).

Capsid DNA extraction and sequencing

Lysates were mixed with PEG8000 1 M NaCL 10% (final concentration) to a final volume of 250 mL and let precipitate overnight at 4 °C. The lysates were then centrifuged at 20,000g for 2–3 h, and the pellet was resuspended by pipetting in 1 mL of DNase buffer and then placed in a shaker at 37 °C for 30 min to fully resuspend the pellet. Then 96 µl high titer purified phage were mixed with 2 µl of TURBO DNAse (Thermo Fisher, REF #AM2238) and 2 µl of benzonase (Abcam, REF #ab270049) and incubated at 37 °C for 1 h. The DNase was then inactivated by adding EDTA 15 mM (final concentration) and a heat-treatment for 10 min at 75 °C. Then, a proteinase K treatment (3,7 mg/mL final concentration) enhanced with SDS (0.74% final concentration) at 55 °C for 45 min was performed to degrade phage capsids. For extraction of the remaining DNA, the PureGene DNA extraction kit (Qiagen) with the protocol for Gram-negative bacteria DNA extraction starting from step 8 (adding directly 100 µl of protein precipitation solution to our 100 µl sample) was used. An overnight DNA precipitation step was performed at 4 °C and the ethanol washes were performed on the next day. DNA was then resuspended for 1 h at 60 °C.

Library preparation and sequencing

For the jmh43 capsids sample, libraries were generated manually following the manufacturer’s protocol for Illumina TruSeq Nano DNA low-throughput library preparation kit (Illumina, San Diego, USA, REF #20015964) with a PCR step of 14 cycles. Briefly, samples were normalized to 500 ng DNA and sheared to 350 pb by sonication with a Covaris E220. AMPure XP beads were used for cleanup and size selection, and then adapters were ligated. Fragment sizes for all libraries were measured using a fragment Analyzer (Agilent) and pooled prior to sequencing. Finally, 1200 pM of the pooled libraries were sequenced on Nextseq2000 Illumina sequencer using a P1 flowcell. Sequencing was performed as pair-end 300 cycles.

For the jmh43::dcas9-repc_HP capsids sample, the Novogene NGS DNA Library Prep Set (Cat No.PT004) was used for library preparation. Index codes were added to each sample. Briefly, the genomic DNA was randomly fragmented to size of 350 bp. DNA fragments were end polished, A-tailed, ligated with adapters, size selected and further amplified by rolling circle amplification. PCR products were then purified (AMPure XP system), their size distribution was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA), and they were quantified using real-time PCR. Libraries were sequenced on an Illumina Novaseq 6000 S4 flowcell with a PE150 strategy.

Capsid DNA analysis

Raw reads from the viral faction (data available at ArrayExpress, accession: E-MTAB-14979) were mapped to the jmh43 chromosome and prophage genome using bwa mem v. 0.7.17-r1188. Samtools v.1.18 was used to obtain bam files and sequence depth [75]. Depth normalization was performed by dividing by the total number of mapped reads on that region.

Capsid concentration by ultracentrifugation

Overnight cultures were diluted to an absorbance of OD₆₀₀ = 0.05 and grown anaerobically in 5 mL cell culture tubes at 37 °C overnight. The next morning, the cultures were centrifuged at 7142g for 15 min at 4 °C, the supernatant was then filtered using 0.45 μm filters. 1.5 mL of each filtrate were pelleted by ultracentrifugation in Beckman coulter 1.5 mL tubes (REF #357448) using an Optima max XP ultracentrifuge with a fixed angle rotor (TLA-55) for 30 min at 50,000 RPM. The pellets were then resuspended in phage buffer (150 mM NaCl, 40 mM Tris-Cl, pH 7.4, 10 mM MgSO₄) and subsequently DNAse treated and heat-treated and used for the ddPCR.

Phylogenetic analyses of Bacteroides and hankyphage-like prophages

Both datasets of 26 hankyphage genomes and 161 Bacteroides genomes were separately processed using the PanACoTA pipeline, version 1.4.0 [31]. The “prepare” step was executed with parameters –min 1e06 (to discard very similar genomes, based on mash distance), –max 0.8 (to discard very different genomes, based on mash distance) and –l90 380 (to discard genomes with high L90 values); the “annotate” step was executed with the parameter –l90 380 (same as before); the “pangenome” step was executed with the parameter -i 0.8 (percentage identity to cluster two proteins in the same family is 80%); the “persistent” step was executed with the parameter -t 0.9 (to select persistent genes found in at least 90% of the genomes); and the “tree” step was executed with the parameter -b 1,000 (to compute the tree using 1,000 bootstraps). The remaining parameters across the different steps of the pipeline were left as default. After the filtering steps of PanACoTA, all hankyphage genomes were kept, and 159 out of the initial 165 Bacteroides genomes were kept. The respective phylogenetic trees were rooted at their midpoint, displayed using iTOL [87] and posteriorly edited to include the tanglegram. Tree files for Bacteroides and hankyphages and their respective IDs can be found in supplementary S17 Data, S18 Data and S4 Data respectively, also available at https://figshare.com/s/3ff18cc2f6cc1edab0ae.

Analysis of metaviromes

Metaviromes with more than 10,000 reads mapped to the reference hankyphage p00 [10] were downloaded and assembled with SPAdes v. v3.15.5 (–metaviral) [88]. Metaviromes containing the hankyphage in less than three contigs were selected for further analysis. Reads from each virome were then mapped to its correponding hankyphage using bwa mem v. 0.7.17-r1188 and the sequence coverage estimated with samtools depth [75].