Murine colitis reveals a disease-associated bacteriophage community

A mouse model of intestinal colitis for the isolation of virus-like particles (VLPs).

To determine if colitis alters the intestinal phage community we studied a mouse model of T cell adoptive transfer that exhibits chronic inflammation of the colon (Fig. S1a)^24,25. Rag1^−/− mice were injected with CD4⁺ CD45RB^High T cells or saline, housed in individual cages and monitored for six weeks (42 days) for weight loss and inflammatory markers (Fig. S1b & S1c). Feces were collected from each mouse at day 0, 14, and 42 (Fig. S1a), and genomic DNA from the whole metagenome and virus-like particles (VLPs) was isolated. Whole metagenome and VLP DNAs were sequenced without amplification. Similar to Kang et al., the resulting viromes are among the first quantitative viromes available from the microbiota²⁶. A similar number of VLPs were recovered from feces when counted by epifluorescence microscopy (Fig. S1d). Bacterial composition analysis of the fecal microbiota showed an increased abundance of Proteobacteria at 42 days (Fig. S1e). Similar results were obtained in mice co-housed with their treatment groups for 30 days after T cell transfer (Fig. S2a). These data support prior reports showing that Proteobacterial communities expand during intestinal inflammation^27–29.

Generation and verification of de novo VLP contigs

To alleviate the limitation of a priori knowledge of virus gene sequence features during virus contig identification we developed a sequence independent approach. Our method employs the quantitative comparison of contigs between whole metagenomes and VLP metagenomes, making use of information obtained from physical VLP enrichments. We used sequence based approaches to validate these contigs as viral (see Materials and Methods). We used a quantitative metagenomic approach to detect viruses whose abundances differ between healthy and inflamed mice using an established method of read recruitment to the de novo assembled contigs (Fig. S3)^26,30. Although we used a read mapping ratio approach to limit contigs of potential bacterial origin within our VLP contig dataset, we cannot rule out the presence of bacterial genomic DNA sequences in viral contigs originating from induced prophages. Our custom VLP database consisted of 1104 contigs ranging in size from 2 kb to 486 kb. Approximately 80% of contigs were 2 to 25 kb in length, 16% were 25 to 100 kb, and the remaining 4% were >100 kb.

To validate the viral origin of the contigs we employed two separate approaches. First, we compared their open reading frames to a set of 25,281 viral protein families (VPFs) from known viruses described in Paez-Espino et al³¹. We identified 278 viral contigs with >70% of genes belonging to VPFs (Table S1). We found 600 and 118 viral contigs with 35–70% and less than 35% of their genes belonging to VPFs, respectively (Table S1). The remaining 108 contigs contained genes that did not match any VPF. Contrary, to the standard application of VPFs for the extraction of viral sequences from metagenomes, which requires meeting stringent criteria, we considered our viral contigs validated if they contained at least one VPF. We further validated our VLP contigs using VirSorter, a tool that identifies viral sequences from complex microbial DNA samples³². VirSorter identified 156 of the VLP contigs as phages and of these, 152 were also identified using VPFs (Table S2).

We next compared our contigs to the IMG/VR viral database, which contains over 264,000 metagenomic viral contigs and 8,000 viral isolate genomes³³. We used a sequence-based classification framework to group closely-related viral sequences³¹. Of the 1104 contigs in our database, 382 viral contigs (representing 34% of the total contigs) grouped within 142 viral clusters (Table S3). All viral contigs from IMG/VR that clustered with the contigs from this study originated from intestinal viruses (Table S3). The remaining 722 viral contigs did not group with any viral contigs in IMG/VR or our curated database and were considered singletons.

Quantitative metagenomics reveals alterations to phage communities during colitis.

To determine whether compositional changes occur within the phage community during colitis, we mapped our VLP reads to all publicly available Caudovirales genomes downloaded from NCBI, all viral isolate and metagenomes in IMG/VR and our curated VLP contig database.

Mapping to phage reference genomes from NCBI showed divergence of the phage community in T cell-treated animals at 42 days, suggesting that disease status correlates with an altered fecal phage community (Fig. 1a). Enterobacteriaceae-specific phages were enriched in T cell-treated mice compared to healthy controls at day 42 (Fig. 1b and Table S4), including phages resembling those that infect Proteobacteria that expand during inflammation (Fig. S1f & S2b). Reads that mapped to phages predicted to infect Enterococci were also more abundant in T cell-treated animals on day 42 (Fig. 1b and Table S4). These data correlated with an increased proportion of phage reads that mapped to known viruses in colitic mice (Fig. S4). The presence of some phages correlated with an increased abundance of prospective host bacteria, suggesting that phage interact with known intestinal pathobionts during colitis. (Fig. S1f & Fig. S2b). In addition, VLP reads that mapped to Caudovirales phage families, including the Siphoviridae, Myoviridae and Podoviridae, were elevated in T cell-treated animals (Fig. 1c). These data accord with the altered phage abundances in human IBD patients²¹. The VLP reads that mapped to the Caudovirales genome dataset accounted for less than 0.01% of the total VLP reads per group of mice, indicating that most of the phage specific VLP reads are unique to our dataset or fall below our identity cutoff.

To improve VLP read mapping for our samples we mapped the VLP reads to the IMG/VR database, as well as the 1104 viral contigs extracted from assemblies of our dataset. IMG/VR contains viral isolate genomes and metagenomes from ocean, soil, sewage and animal hosts³³. A total of 1068 viral genomes from IMG/VR recruited reads from our mouse VLP samples (Table S5). 958 of these viral genomes originate from intestinal samples including 406 from humans, 417 from mouse, 102 from rat and 33 from other mammals (Fig. S5 & Table S5). This increased our read mapping percentage by two orders of magnitude (1.2%).

We next mapped each sample’s VLP reads to our curated database of 1104 viral contigs (Fig 2a). Using the curated contig dataset allowed us increase our read mapping average to 22% (Table S6), yet 78% of the reads were not mapped to any database. Comparative read mapping analysis of day 0 samples showed both T cell-treated and control animals had a similar viral community. However, as colitis progressed the viral community shifted, with the most marked shift occurring 42 days after T cell treatment (Fig. 2a, 2b & 2c). This shift in viral community composition was also observed in animals that remained co-housed with their treatment group for the duration of the experiment (Fig. S6). Thus, T cell-mediated colitis was associated with a divergence of the intestinal viral community.

The shift in phage community composition at day 42 reflected a decrease in viral contig abundances (Table S7) and a significant loss of diversity leading to an increased dissimilarity between healthy and diseased mice (Fig. 2c). Rarefaction analysis revealed sufficient read mapping coverage for all samples, ruling out low sampling as a cause of the divergence (Fig. S7). In addition to altered phage abundances there was a marked divergence in the bacterial community during colitis (Fig. S8). We found that some viral contigs differed significantly between control and diseased mice at day 42 (Fig. 2d). Of the contigs that were statistically significant, 35 contigs were at least 10-fold more abundant and 41 contigs were at least 8-fold less abundant during colitis.

When comparing the presence of various contigs between healthy and colitic animals, we observed a stochastic pattern consistent with a random distribution of contigs within animals at various time points (Fig. 3a). 4% of the VLP contigs were shared between healthy animals as compared to only 1% in animals with colitis (Fig. 3a). Furthermore, the colitic animals contained a greater number of singleton viral contigs (Fig 3b). Although contigs shared between healthy animals accounted for a small proportion of the total VLPs, these data suggest that the viral community diverges as the intestine transitions from a healthy to dysbiotic state.

At day 42 there was a dramatic reduction in the number of reads recruited to VLP contigs in colitic mice. This was most pronounced for mouse number two (Fig. 2a). 20% of the mapped reads from this animal aligned with four highly abundant contigs. Three of these contigs were among the most differentially abundant contigs at the height of disease (Fig. 2d). An 18.4 kb contig recruited 8.3% of the reads and the other three contigs of 136.6 kb, 145.5 kb and 176.1 kb recruited the remaining 11.7 % of the reads (Fig. S9 & Table S8). Although these contigs recruited a large proportion of the VLP reads from T cell-treated animals, read recruitment to these contigs varied between animals (Fig. S9).

We were able to assign phage taxonomy or bacterial host information to two of these four contigs. We assigned phage taxonomy or host information using a list of signature phage orthologous groups (POGs), phage-encoded tRNAs, and clustered regularly interspaced short palindromic repeat (CRISPR) spacers that serve as a historical record of phage-bacteria interactions^34–36. The 145.5 kb contig contained a tRNA 100% identical to the genus Acholeplasma, a member of the Mollicutes. This contig had a region of ~8 kb that resembled known phage genes (Fig. S9 & Table S8). The 136.7 kb contig matched a POG belonging to Spounaviridae, a subfamily of the Myoviridae composed of phages known to infect pathogens such as Staphylococcus aureus³⁷.

Of the remaining 35 statistically significant contigs of high abundance, only three other contigs could be assigned taxonomy or host information (Table S9). A tRNA broadly connected one contig to either the Firmicutes or Tenericutes, one contig was assigned to a POG representing the genus Betalipothrixviruses which infect Archaea³⁸, and the third contig was linked to the Clostridiales through a CRISPR spacer. Although we assigned taxonomy or host information to only a few significantly different contigs, this was expected based on previous studies^21,39. This highlights the power of our approach for identifying VLPs independent of existing database information. Of the 41 contigs with significantly reduced abundance, CRISPR spacer and tRNA matching revealed six to be potential Clostridiales phages and two were identified as phages that infect the Lachnospiraceae family within the Clostridiales. These data indicate that colitis is associated with a reduction in Clostridiales-like phages.

VLP identification reveals changes in the phage community during colitis

We next determined taxonomic or host assignments for all 1104 contigs in the curated VLP database. We were able to assign 116 contigs (11%) to a specific phage taxonomic group or host designation using one or more of the three identification methods described above. The number of total assigned contigs is smaller than the number of contigs predicted to be of viral origin using VPFs because numerous VPFs only contain viral hallmark genes and lack taxonomic or host information. CRISPR spacer matching resulted in the greatest number of assigned contigs (50 contigs), followed by tRNAs (41 contigs) and POGs (36 contigs) (Fig. S10). Both tRNA and CRISPR spacer matching identified 11 shared contigs, accounting for 27% and 22% of the total contigs identified using those methods respectively. Of the 36 contigs identified using POGs, 26% belonged to the Myoviridae, 18% to the Siphoviridae and 37% to the Podoviridae. The remaining 19% of contigs matched several archaeal viruses belonging to the Fuselloviridae and Lipothrixviridae (Fig. S10a).

POG matching identified numerous phages with taxonomic signatures of known Enterobacteriaceae-like phages. This is consistent with the increased abundance of Proteobacteria in the metagenomes of T cell-treated mice and the mouse VLP reads mapping with higher abundance to Proteobacterial phage genomes (Fig. S1e & Fig. 1a). Since viruses within a given POG can infect very different bacterial hosts, we cannot completely exclude that the contigs classified as Enterobacteriaceae-like phages are from viruses that infect bacteria outside the Proteobacteria. Additionally, POGs identified phages taxonomically related to phiCD119, phiKZ, phiHp1, AHJD and Spounavirinae phages that infect pathobionts including Clostridium difficile, Pseudomonas aeruginosa, Haemophilus influenza, and Staphylococcus aureus^37,40–42 (Fig. S10a & Table S9). We were able to obtain order, family and in some cases genus-specific identification of multiple phage contigs using tRNA and CRIPSR spacer matching. As expected, the most highly represented phages in our database matched hosts in the phyla Bacteroidetes and Firmicutes. Firmicute-specific phages were dominated by those predicted to infect Clostridiales bacteria including phages specific for Lachnospiraceae, Ruminococcaceae and Clostridia. The next most prevalent were phages specific for the Streptococci. Bacteroidetes phages were almost exclusively represented by Alistipes phages (Fig. S10b, S10c & Table S9).

To identify differences in the abundance of host-assigned phage contigs we grouped the contigs based on host identity, summed the read-mapping derived abundances and averaged these values for healthy and colitic animals. At the order level there was a significant decrease in the abundance of Clostridiales phages in T cell-treated animals (Fig. 4a). This supports our finding, from the analysis of the differentially abundant contigs that a2 reduction in Clostridiales phages occurs in colitis. When these data were stratified into family- and genus-specific contigs we discovered that phages predicted to infect Clostridium sp. KLE 1755 and some Lachnospiraceae bacteria decreased, whereas Lachnospiraceae phages that target strain COE1 were more abundant during colitis (Fig. 4b, 4c & 4d). Additionally, Clostridiales phages predicted to infect Clostridium difficile were more abundant in the T cell-treated animals (Fig. 4e).

Streptococcus sp. phages were more abundant in T cell-treated animals, supporting a prior report on intestinal phages from IBD patients²¹, as well as phages predicted to infect Alistipes (Fig. 4f & 4g). The 16S rRNA gene read abundances for the host bacteria infected by these phages were not statistically different between healthy and colitic animals, although Alistipes trended toward lower abundance during colitis (Fig. 4h). This suggests that shifts in phage abundance are not solely dependent on bacterial host density. Together, these data support the idea of a defined phage community during colitis.

Mouse VLP reads map to phages from human phageomes

To determine if mouse VLP sequencing reads matched intestinal phages from humans, we mapped our VLP reads from healthy and colitic animals to groups of phage contigs derived from healthy humans and humans with IBD^21,39. A total of 38 contigs, accounting for 25% of the healthy human phageome, recruited reads from both healthy and colitic animals (Fig. 5a). Comparative analysis of the read mapping frequency revealed multiple healthy human phage contigs that recruited fewer reads from colitic animals and seven phages that were exclusively in healthy animals (Fig. 5a). Consistent with the observation that more contig connections were maintained under healthy conditions across time points (Fig. 3a), these data suggest that the healthy component of the murine phageome overlaps with that of humans and that the abundances of these phages decrease during colitis.

We next assessed the abundance of VLP reads that mapped to a set of 519 VLP contigs generated from the Chicago cohort of Crohn’s and ulcerative colitis patients described by Norman et al.²¹. A total of 169 contigs, accounting for 33% of VLP contigs generated from IBD patients, recruited reads from either healthy or colitic animals (Fig. 5b). Phages from diseased humans recruited more reads from colitic animals and in some cases only human IBD phage contig specific reads were represented in colitic mice (Fig. 5b). Thus, the murine phageome overlaps with that of humans both during intestinal health and disease.

Animals.

Two-week-old male C57BL6/J Rag1^−/− mice⁵¹ were weaned and co-housed together in the same cage for four weeks. At six weeks of age three mice C57BL6/J Rag1^−/− mice were administered 5×10⁵ naïve CD4⁺ CD45RB^High T cells via retro-orbital injection. A second group of three C57BL6/J Rag1^−/− mice were given a saline injection as a control. All mice were housed in either separate cages or by treatment group for the duration of the experiment. The sample size of three mice per group was used because this was the minimum number of mice needed for reliable statistical analyses. Four male wild-type C57BL6/J mice aged six to eight weeks were used as spleen donors for T cell isolation. Male germ free C57BL6/J Rag1^−/− mice were reared in sterile isolators as described previously⁵². Fecal samples from three eight week old germ free Rag1^−/− were pooled in duplicate. After co-housing for four weeks experimental mice were randomly assigned to treatment groups. It was not possible to blind researchers to the treatment conditions of the animals prior to experimentation. The Institutional Animal Care and Use Committee of UT Southwestern Medical Center approved all animal protocols.

T cell isolation.

The spleens from four male wild-type C57BL6/J mice were dissected and added to RPMI (BD Biosciences, San Jose, CA) containing 10% fetal bovine serum (Gemini Bio-Products, Sacramento, CA). The spleens were homogenized using the frosted ends of sterile microscope slides and passed through a 70 μm cell strainer. CD4 cells were negatively selected using a Dynabeads Untouched Mouse CD4 Cells Kit (Invitrogen, Carlsbad, CA). Purified mouse CD4 cells were Fc-blocked using purified rat anti-mouse CD16/CD32, clone 2.4G2, (BD Biosciences) followed by two color staining using anti-mouse CD4-FITC, clone RM4–5 and anti-mouse CD45 RB-APC, clone C363–16A (Thermo Fisher, Waltham, MA). CD4⁺ CD45RB^High cells were collected on a MoFlo XDP cell-sorting machine at the UT Southwestern Flow Cytometry Core Facility.

Metagenomic DNA isolation.

VLPs and their DNA were isolated as described previously⁵³. Briefly, fecal pellets were homogenized in 2 ml of SM+ buffer⁵⁰ and 1.3 ml sample volume was clarified by centrifugation (2500g for 5 min followed by 5000g for 15 min at 4ºC). The supernatant was filtered through a 0.45 ¼m Millex-HV syringe filter (EMD Millipore, Burlington, MA) and treated with 100 units of DNase I (Roche, Pleasanton, CA) and 15 units of RNase A (Sigma, St. Louis, MO) for 1 hr at 37ºC. VLPs were not further purified by cesium chloride (CsCl) gradient centrifugation. Although CsCl purification is efficient at removing potential contaminating nucleic acids, this method is not suited for quantitative metagenomic comparison of a large number of samples due to its limited throughput⁵³. Additionally, this method discriminates against specific phages and can vary in reproducibility due to the necessity of removing the VLPs from precise locations within the gradient that are difficult to normalize between samples⁵³. Therefore, we proceeded to isolate VLP DNA from filtered and DNase/RNase treated samples by adding 100 mM EDTA, 50 ¼g/ml Proteinase K (New England Biolabs, Ipswich, MA) and 0.5% sodium dodecyl sulfate to each sample and to incubate them at 56ºC for 30 min. The samples were extracted with an equal volume of phenol:chloroform:isoamyl alcohol 25:24:1 (Sigma) followed by an additional extraction using just chloroform. The DNA was precipitated with isopropanol and cleaned using a MinElute Reaction Cleanup Kit (Qiagen, Hilden, Germany). Following a rigorous VLP purification strategy can greatly reduce the possibility of contaminating prokaryotic host genomic DNA that can lead to false positive viral gene identification from the VLP fraction^54,55. For the isolation of total metagenomic DNA, the remaining 0.7 ml of fecal material was treated with 5 mg/ml lysozyme (Amersco, Solon, OH) at 37ºC for 30 min. The samples were transferred to 2 ml Lysing Matrix B tubes (MP Biomedicals, Solon, OH) and were disrupted by bead-beating in a Bullet Blender 5 (Next Advance, Troy, NY) at a speed setting of 9 for 1 minute followed by a 2 minute incubation on ice. This was repeated 5 times. Samples were centrifuged at 16,300g for 1 minute and DNA from the supernatant was extracted and purified as described above.

Metagenomic DNA sequencing.

DNA sequencing was performed using Illumina HiSeq 2500 Rapid Run at 150 cycles in paired-end mode. DNA concentration and quality was determined using the Qubit dsDNA HS kit (Invitrogen) and a Bioanalyzer 2100 (Agilent, Santa Clara, CA). DNA libraries were prepared using a KAPA High Throughput Library Preparation Kit with standard PCR library amplification (KAPA Biosystems, Wilmington, MA). TruSeq adapters (Illumina, San Diego, CA) with 6 bp indices were used to enable multiplexed sequencing. Illumina libraries were size selected using AMPure XP beads (Beckman Coulter, Indianapolis, IN). Library quality was measured on a Bioanalyzer 2100 by product size and concentration. On average, 43 million paired-end whole metagenome reads and 30 million paired-end VLP reads were generated for each sample. Sequencing reads were demultiplexed using the 6 bp TruSeq adapter sequences allowing for one mismatch using the Illumina CASAVA software. DNA quality control, HiSeq 2500 sequencing and library preparation was performed by the UT Southwestern Eugene McDermott Center for Human Growth and Development Next Generation Sequencing Core.

Decontamination and processing of reads.

Read decontamination and trimming was performed using the BBMap short read aligner v36.99 (https://sourceforge.net/projects/bbmap). Raw reads were mapped to the internal Illumina control phiX174 (J02482.1) and the mouse (mm10) and human (hg38) reference genomes using the bbsplit algorithm with default settings. The resulting unmapped reads were adapter trimmed and low-quality reads and reads of insufficient length were removed using the bbduk algorithm with the following parameters: ktrim=lr, k=20, mink=4, minlength=20, qtrim=f. These reads went through a second round of decontamination by mapping to a set of reference contigs assembled from samples obtained from the feces of germ-free C57BL6/J Rag1^−/− mice using bbsplit. These final unmapped reads constituted the cleaned M and P read sets.

Read based quantification of taxa.

To assess microbial community composition we used phyloFlash Version 2.0 (https://github.com/HRGV/phyloFlash) to reconstruct 16S and 18S rRNA genes from the Illumina read datasets and quantified taxa by read mapping. 16S and 18S rRNA genes were classified with phyloFlash against the SILVA SSU Ref NR 99 database 123⁵⁶. Taxa abundances were summed at the phylum level and the genus level for specific genera.

Metagenomic assembly.

Whole metagenome, VLP and germ free mouse decontamination reference assemblies were performed using the IDBA-UD assembler v1.1.1⁵⁷. Prior to assembly, paired end reads were interleaved using the fq2fa merge filter command within IDBA-UD. Assemblies were generated in the absence of error correction using the default parameters.

Construction of the curated VLP contig database.

VLP reads were assembled into 18 sets of contigs, with one set per animal per time point. VLP contigs smaller than 2 kb were removed from each sample’s contig set. Using these pruned contig sets, each sample’s VLP and whole metagenome reads were mapped to their corresponding sample’s contigs using the bbmap algorithm within BBMap with the following settings: ambiguous=random, qtrim=lr, minid=0.97. Read mapping coverage for each contig was generated by setting the ‘covstats’ parameter. The percentage of reads that mapped to each contig was used to calculate a read-mapping ratio per contig (VLP/whole metagenome). We assumed that if this ratio was greater than 1 (i.e., more VLP reads than whole metagenome reads mapped to a contig), then these contigs were more highly represented in the viral fraction and likely represented genuine viral contigs. From each of the 18 read mapping sets we selected and pooled the top 100 contigs based on the largest ratios. We removed redundancy by sequence clustering using cd-hit-est⁵⁸ at an identity threshold of 99%, resulting in a collection of 1104 unique contigs that constituted our curated viral contig database.

Quantitative read mapping.

The complete genomes of 1993 Caudovirales genomes downloaded from NCBI, a group of 154 phage genomes found in healthy humans from the Manrique study³⁹, contigs obtained from a Chicago cohort of IBD patients from the Norman study²¹ and our curated viral contig database generated from healthy and colitic mice were used as mapping references using the bbmap algorithm. For the Caudovirales genomes and the curated viral contig database we used the following bbmap parameters: ambiguous=random, qtrim=lr, minid=0.97. For the Manrique and Norman datasets we decreased the identity threshold to minid=0.90, to account for increased heterogeneity between mouse and human samples³⁹. Read mapping statistics for each contig were generated by setting the ‘scafstats’ parameter. For the Caudovirales genomes, the Manrique dataset and the Norman dataset the total number of unambiguous reads per contig were tallied and averaged across animal treatments and time points. For the curated viral contig database, the percent ambiguous and unambiguous mapped reads were combined for each contig and compared between treatments and time points.

Confirming viral contigs, establishing level of similarity with existing viruses, and defining host associations.

We used a curated set of 25,281 viral protein families (VPFs)³¹ and the VirSorter database³² to compare open reading frames within the 1104 viral contigs from this study for viral protein similarity. VPFs are a recently developed tool yet to be widely used for virus discovery. VPFs have been benchmarked for high specificity for viruses and only in rare cases a VPF may include a prokaryotic host protein. Additionally, we obtained taxonomic information for viral contigs by using phage orthologous gene clusters (POGs) that are recognized as taxon-specific signatures for virus classification³⁵. For all Hidden Markov Models searches we used the algorithm hmmsearch with trusted cutoffs from the hmmer web server⁵⁹ and an e-value below 1e^-05. To cluster viral sequences we used the public viral database IMG/VR composed of over 264 thousand metagenomic viral contigs and almost 8 thousand virus isolates³³. We used the sequence-based classification framework described in Paez-Espino et al.³¹ for grouping closely related viral sequences. The method relies on both nucleotide identity and the fraction of total alignment for pairwise comparisons of viral sequences. We used BLASTn in the Blast+ package⁶⁰ to find hits with an e-value cutoff of 1e⁻¹⁰ and at least 50% identity across ≥ 50% of the shortest sequence length. We used three different computational methods to predict host-virus connections: i) Grouping with isolated viruses, metagenomic viral contigs and host or taxonomic association using POGs. ii) Prokaryotic CRISPR spacer matching to viral contigs requiring at least 95% identity over the whole spacer length. To acquire CRISPR spacers from our experimental dataset each animal’s whole metagenome reads from each time point were combined and contigs were assembled using IDBA-UD. Binning of metagenome assemblies to retrieve genome bins was performed using MetaWatt⁶¹ with default settings except avg read length=75, contig length threshold=500 and minimum bin seed size=5000. Taxonomy was assigned to each bin based on MetaWatt taxonomy binning. CRISPR repeat patterns and spacers were extracted from bins using metaCRT⁶² with the following settings minNR=4, minRL=19, maxRL=48. Spacers were extracted from CRISPRs and aligned to the curated viral contig database using the seal algorithm in the BBMap software suite with the k parameter set to 19. iii) Matching viral transfer RNA sequences to bacterial transfer RNAs. We identified tRNAs from our viral contig set using ARAGORN v1.2⁶³ and compared these against isolate microbial genomes in IMG/M⁶⁴ using BLASTn. Only hits with a 100% identity over 100% of the length were considered. There are multiple methods available for phage host identification from VLP metagenomes⁶⁵. We chose to use CRISPRs and tRNAs to make host connections because CRISPR matching does not depend on preexisting databases for sequence comparison and tRNAs are common features of phage genomes that provide host specific information without the need for a reference phage genome.