Assembly and initial analysis

The average read length in our VLP data set was longer than in most previous metagenomic studies of viral DNA, allowing extensive assembly of individual reads into contigs. We assembled VLP sequences using the Newbler assembler (Miller et al. 2010) (40 bp overlap, 90% identity), which yielded 7175 contigs at least 500 bp in length. Fully 86.6% of the sequence reads were recruited into these contigs (Fig. 2A). The longest contig was 46 kb, 73 contigs were longer than 10 kb, 279 contigs were longer than 5 kb, and 3028 contigs were longer than 1 kb.

Open in a separate window

Figure 2.

Assembly and functional annotation of shotgun metagenomic sequences from the human gut virome. (A) Analysis of recruitment of VLP sequence reads into contigs. The y-axis shows the number of sequence reads, the x-axis shows contig length. Pyrosequencing data from the human virome were assembled into 7147 contigs up to 47.8 kb in length. Linear contigs are shown in yellow, circular contigs are shown in blue. (B) Analysis of protein functions in VLP contigs. The functions encoded in VLP contigs were predicted using the Pfam database, then grouped using custom database-relating Pfam domain identifiers to phage functions (Supplemental Table S4). The relative proportions of pyrosequencing reads falling within ORFs of different annotations were plotted according to their sample of origin (y-axis). Each bar is indicated by the sample code, where L or H indicates low- or high-fat diet, the adjacent number indicates the subject number, and the number after the hyphen indicates the day of the study. “X-1” indicates ad-lib diet. (C) Taxonomic classification of VLP communities is consistent across samples. Samples (in columns, labeled as in Figs. 1, ​,5)5) are characterized according to the number of sequences from each sample that are assembled into a contig that is classified by taxonomic family. Phage families are indicated by the color code to the right. “Unknown” (black) indicates contigs that cannot be classified in any way. “Multiple hits” (white) indicate contigs that have proteins that are similar to multiple families.

The approximate size of the VLP community can be estimated by PHACCS (PHAge Communities from Contig Spectrum) (Angly et al. 2005), which calculates the degree to which a group of sequences coassemble into contigs, and compares them with simulated communities of different sizes. Similar to what has been seen recently in the human gut virome (Reyes et al. 2010), the median richness (number of species) of the 16 samples was 44 (range = 19–785), and the average Shannon Diversity was 3.46 (SD = 0.59). The most abundant genotype was predicted to account for 16.20% of the total (SD = 2.09%), predicting the complete assembly of the most abundant VLP genomes in our study.

Analysis of gene content

To characterize these VLP contigs, we compared both nucleotide sequences and open reading frames [ORFs] with (1) the NCBI nonredundant database, (2) the Pfam database of conserved amino acid motifs (Sonnhammer et al. 1998), (3) the Clusters of Orthologous Groups (COG) database of annotated bacterial protein families (Tatusov et al. 2003), (4) A CLAssification of Mobile genetic Elements (ACLAME) (Leplae et al. 2009), (5) the Antibiotic Resistance Genes Database (ARDB) (Liu and Pop 2009), and (6) the Virulence Factors Database (VFDB) (Yang et al. 2008). All VLP contigs and associated annotation can be viewed using a web-based interface at http://microb215.med.upenn.edu/cgi-bin/gb2/gbrowse/phage_metagenomics/. A total of 22% of ORFs contained recognizable Pfam motifs. Essential bacteriophage functions were well represented, including functions required for both lytic and lysogenic growth (Fig. 2B).

VLP contigs were classified according to their similarity to ICTV-defined bacteriophage families (Fig. 2C). There were 1268 contigs with amino acid similarity to members of the Siphoviridae family (18% of the total), 686 (10%) to Myoviridae, 344 (4.8%) to Podoviridae, 68 (0.9%) to Microviridae, and 0.4% to other families. Of the remaining contigs, 813 (11%) had amino acid similarity to multiple bacteriophage families, and 3969 (55%) did not have significant similarity to any bacteriophage families. Membership in these families was similar among the subjects studied (Fig. 2C). No strong candidates for eukaryotic viruses were detected either through nucleotide comparison to known eukaryotic viral genome sequences, or through similarity to conserved eukaryotic protein domains. In a few cases, a gene was annotated as similar to a eukaryotic virus, but for each of these genes contigs could be identified that also encoded multiple phage genes, suggesting that these proteins contain motifs common to both bacterial and eukaryotic viruses. Thus, this classification indicates that the viruses studied here were comprised primarily of tailed bacteriophage, with some representation of additional DNA phage families.

Analysis showed that nine VLP contigs formed closed circles, 4.5–6 kb in length, suggestive of completion of these genome sequences. None of the nine had significant nucleotide similarity to any sequence in the NCBI nonredundant database (E-value < 10⁻³). However, eight of the nine contigs contain an ORF that aligns to the capsid F protein sequence from Microphages (Rohwer and Edwards 2002), the family of 4.5–6-kb circular ssDNA phage containing the prototype bacteriophage X174. Three of the eight contigs were closely related to the chp1-like Microphage. Those eight genomes also contain proteins similar to phage proteins involved in replication (n = 4), proteolysis (n = 2), and scaffolding assembly (n = 2). Only two of these genomes have significant sequence similarity to each other—contigs c03390 and c04421 are 93% identical, but were too diverse to coassemble. The ninth genome (contig c04570) contains proteins that resemble those found on plasmids from a wide range of bacteria, primarily Firmicutes, and may represent a temperate phage that maintains itself as a plasmid (Sternberg and Austin 1981).

Comparison of total community and VLP metagenomes

We next investigated the proportion of phage DNA in our total (bacterial) metagenomic data set. We calculated the proportion of VLP sequences that have a significant match in the total shotgun data set (E-value < 10⁻⁵), and vice versa (Fig. 3A). As expected (Reyes et al. 2010), the VLP sequences represent a minority of the total DNA from stool, in the range of from 4% to 17%, although this value is dependent on sampling depth, because the VLP community is a subset of the total community. We also assembled and annotated the total community shotgun sequence data set (347 Mb) using COG (Tatusov et al. 2003), then compared it with the VLP data set. The proportions of COG classes present in the total community (mostly bacteria) and VLPs were quite different, as expected (Fig. 3B). Bacteria were significantly enriched in genes for synthesis of amino acid and carbohydrate precursors, ion transport and metabolism, translational machinery, and cell wall/membrane biogenesis, while VLP contigs were enriched in genes for replication, recombination and repair, and unknown functions (P < 0.05, t-test). Thus, the profile indicates that viruses recruit host cell machinery for translation, energy production, and synthesis of macromolecular precursors while using their own coding capacity to encode functions for replication.

Open in a separate window

Figure 3.

Comparison of gene content in total microbial communities and VLP communities. (A) The proportions of total genes identified in shotgun metagenomic analysis of total stool DNA (left) is compared with genes in VLP communities (right). (B) VLP DNA encodes biochemical functionality that is markedly different from the total microbiome. Function annotation was performed according to comparison to the COG database (Tatusov et al. 2003). The proportion of reads from each assembled data set that fall within an ORF of the indicated annotation are plotted on the x-axis (mean ± SE).

Prophage abundance

A key question in investigating the gut virome centers on what fraction of phages are temperate, since this group can install new genes in bacteria and alter phenotype via lysogenic conversion (Brussow et al. 2004; Reyes et al. 2010). We used three indications of a temperate lifestyle to annotate the VLP contigs: (1) nucleotide identity to sequenced bacterial genomes, which is indicative of prophage formation (90% of bases aligned at 90% identity), (2) presence of integrase genes (according to Pfam) (Sonnhammer et al. 1998) and COG (Tatusov et al. 2003) annotations, and (3) significant similarity of multiple proteins to prophages annotated in the ACLAME (Leplae et al. 2009) database of mobile genetic elements (E-value < 10⁻⁵). This strategy provides a minimal estimate of the number of temperate phages in our data set, since authentic temperate phages may not be positive by any of these criteria. Of the 3029 VLP contigs of at least 1 kb in length, 428 (14%) had multiple proteins that significantly resembled an annotated prophage, 73 (2.4%) contained an integrase gene, and 37 (1.2%) aligned to a sequenced bacterial genome (Fig. 4A). Of the 505 contigs that fit at least one of these three criteria, 442 (88%) also contained an annotated bacteriophage gene. When individual VLP reads were mapped to known bacterial genomes, only 1.5% (13,808) aligned at 90% identity or higher. Of those reads, 74% map to just five genomes, Bacteroides vulgatus ATCC 8482 (n = 4555), Eubacterium eligens ATCC 27750 (n = 2147), Faecalibacterium prausnitzii L2/6 (n = 2081), Parabacteroides distasonis ATCC 8503 (n = 785), and Bacteroides thetaiotaomicron VPI-5482 (n = 583). These reads aligned to short, mostly contiguous regions of each genome and assembled into contigs that encode bacteriophage proteins (Fig. 4B,C), as expected for prophages. All of the reads matching Eubacterium eligens (n = 2147) mapped to its plasmid, which also contains annotated phage proteins, suggesting that this plasmid may contain an integrated prophage or itself be a nonintegrating prophage, potentially analogous to phage P1 (Sternberg and Austin 1981). Of our 3029 VLP contigs over 1 kb, a minimum of 505 (or 17%), can be tentatively classified as temperate phages by at least one of the above criteria.

Open in a separate window

Figure 4.

Analysis of temperate bacteriophages in the human gut virome. (A) Venn diagram indicating frequency of functions associated with temperate phages. VLP contigs were annotated according to the presence of integrase-like sequences (orange), BLAST alignment to sequenced bacterial genomes, suggestive of prophage formation (90% length at 90% identity; green), and presence of multiple genes with significant similarity to a prophage element within the ACLAME database (red). Map of VLP pyrosequencing reads aligning to the genomes of (B) Faecalibacterium prausnitzii L2/6 and (C) Parabacteroides distasonis ATCC 8503 (90% length at 90% identity). (Inset) Contigs that correspond to each peak, the reads that make up each contig, and their annotated genes. Top strand reads are indicated by blue, bottom strand reads by gray. (D) Prevalence of prophage sequences according to bacterial phylum of origin. Contigs with amino acid similarity (E-value < 10⁻⁵) to prophage sequences within the ACLAME databases are shown according to which bacterial phylum they match. Any contig with similarity to more than one phylum is classified as “Multiple.” Roughly 39% of VLP contigs (n = 2814) did not have significant similarity to any prophage sequence in this database.

Variation among animals in prophage induction has been previously shown in mouse models (Reyes et al. 2010), leading us to investigate induction here. We probed the level of induction of the above five prophages indirectly by comparing shotgun sequences from total community DNA (mostly bacterial), with VLP sequences (Supplemental Fig. S2). In only one case (F. prausnitzii) did the abundance of the bacterial host correlate with the abundance of its corresponding prophage. The lack of correlation of the other four prophages suggests the possibility that the degree of induction varied among the individuals studied.

We were able to infer the groups of bacteria infected by these putative temperate phage by comparing the VLP contigs to prophage sequences in the ACLAME database. Of the 2814 contigs that had significant amino acid similarity (E-value < 10⁻⁵) to an ACLAME prophage sequence, the contigs that were assigned exclusively to the Firmicutes accounted for fully 41% of the total (n = 1148) (Fig. 4D). Surprisingly, similarity to prophage from the other common gut resident, Bacteroidetes, accounted for <0.9% of the assigned contigs (n = 24). While Proteobacteria account for only a modest fraction of the total gut bacteria, contigs with similarity to Proteobacteria accounted for 16% (n = 447) of the total. This pattern also extended to contigs that were similar to multiple bacterial phyla (“Multiple” in Fig. 4D). Of the 1176 VLP contigs that were similar to more than one bacterial Phylum, 34% had a match to Firmicutes, 0.8% had a match to Bacteroides, and 13% had a match to Proteobacteria. This suggests that the prevalence of temperate bacteriophage in the gut may differ among bacterial phyla, though the conclusions depend on representation of phage sequences in the databases used for analysis, and so conclusions may evolve as further sequence data accumulates.

Functionality in the gut virome: CRISPRs

Clustered Regularly Interspaced Short Palindromic Repeats (CRISPRs) are mediators of a recently described form of bacterial adaptive immunity. Short DNA sequences (26–72 nt in length) are captured from invading DNA elements and installed as CRISPR spacer sequences in bacterial chromosomes. Multiple spacer sequences are separated by partially palindromic repeats within CRISPR arrays. After transcription of CRISPR arrays, RNA copies of the CRISPR spacer sequences act as recognition modules in nucleoprotein complexes, which bind to targeted nucleic acid from genomic invaders such as phage or plasmids, and program their degradation (Brouns et al. 2008; Sorek et al. 2008). Thus, CRISPR spacers provide a record of the genomic parasites that bacteria have encountered.

Analysis of CRISPRs in our assembled data provided a detailed record of phage–host and phage–phage competition. A total of 38 CRISPR arrays were found in shotgun metagenomic assemblies for our total community data set. The 38 contigs containing these CRISPRs aligned primarily to genomes of gut Bacteriodetes (n = 16, E-value < 10⁻⁵). Of the 45 spacer sequences found within those 38 CRISPR arrays, fully 80% (n = 36) showed no matches to any known sequence, but one spacer perfectly matched a VLP contig sequence (contig c05189; 100% identity over 43 bases). Both the bacterially encoded CRISPR spacer and VLP target were found within a single individual, suggesting that the CRISPR spacer sequence may have restricted phage replication within that subject.

Unexpectedly, we also found 22 CRISPR arrays in our VLP contigs. CRISPRs have not previously been reported from free phage particles, but one report did identify CRISPR sequences in potential prophages of Clostridium difficile (Sebaihia 2006). A related CRISPR array was found in our virome samples (with >95% identity in repeat regions to the annotated C. difficile CRISPR5). Analysis of novel CRISPR spacer sequences showed one high-stringency match between a VLP spacer (on contig c05834) and a separate VLP contig (c02690) (95% identity over 39 bases). The CRISPR spacer–target pair were isolated from the same individual and both were detected at the same two time points. Phage are well known to encode an extensive set of functions for competing with other phage (Calendar 1988; Refardt 2011)—these data indicate that CRISPRs may mediate phage–phage competition as well. It will be valuable to assess CRISPR evolution and targeting in further phage metagenomic data sets.

Functionality in the gut virome: Antibiotic resistance

To interrogate antibiotic resistance in the VLP metagenomic data set, the ORFs found on contigs and unassembled reads were compared with known antibiotic resistance genes (BLASTp against the ARDB database) (Liu and Pop 2009) E-value < 10⁻⁵, yielding 614 matches to antibiotic resistance genes. These included multidrug efflux transporters (n = 355), vancomycin resistance genes (n = 129), tetracycline resistance genes (n = 18), and beta-lactamases (n = 16). The 33 highest quality matches to assembled contigs were examined further (Supplemental Table S2). Some of these encoded clear antibiotic resistance genes (e.g., beta-lactamase, drug efflux pumps, streptogramin acetyltransferase), while others encoded relatives of antibiotic resistance proteins that are of uncertain importance to resistance (e.g., ABC transporters and the VanRS two-component signaling system). Five of the contigs also contained identifiable phage genes. Transmissible antibiotic resistance is believed to involve primarily mobilization by plasmids and transposons (Barlow 2009; Juhas et al. 2009)—our results indicate that the role of mobilization by phage may deserve further investigation.

Variation of the gut virome among individuals and during dietary intervention

We next asked how viromes varied among individuals and how the dietary intervention affected virome structure. We characterized each virome sample by enumerating the proportion of sequence reads recruited into each VLP contig over the full contig data set. In Figure 5A the proportional abundance of each VLP contig is shown by the color code. The Euclidean distance was then computed between all pairs of virome samples and used for statistical analysis. Comparison of the collection of within-subject distances among time points to between-subject distances showed that between-subject distances were significantly greater (Fig. 5B; P < 0.001, significance assessed by permutation of sample labels). Thus, each individual contained a unique virome that was globally stable over the 8-d time course. A previous study showed individual distinctiveness and stability of the human virome over a year (Reyes et al. 2010).

Open in a separate window

Figure 5.

Alterations in VLP contig abundance associated with diet. (A) Proportions of VLP sequence reads in contigs from different subjects and time points. Vertical bars indicate the proportion of reads within each contig. The contigs are shown in columns, and subject/time-point combinations in rows. Hierarchical clustering was performed on both rows and columns according to Euclidean distance and complete distance agglomeration. Samples are labeled by subject according to diet: high-fat (H1, H2), low-fat (L1, L2, L3), and ad-lib (X), as well as day of dietary intervention (days 1, 2, 7, or 8). The proportion of all VLP reads in each contig are shown by the scale at the bottom. (B,C,D) The Euclidean distance between samples is shown according to median (line), quartile (box), and range (whisker). (B) Between-subject variation is significantly greater than within-subject variation (P < 10⁻⁴, subject label permutation). (C) The distance between the gut viromes of individuals on the same diet (left) was significantly smaller at the end of their dietary treatment than it was at the start (P = 0.05, permutation of diet labels), while there was no increase in similarity for individuals on different diets (right). (D) The distances within subjects that was measured between samples taken on the first 2 d of the timecourse (left) and the last 2 d of the timecourse (right) were significantly lower than the distances between those two sets of days (middle; P = 0.04, permutation of day labels). (E) Covariation of bacterial and VLP community diversity. Distances between pairs of bacterial and VLP communities were calculated as described in the Methods. The similarity of bacterial and VLP communities is shown through PCoA analysis, where each data set was rotated and scaled for maximum superimposition. Each circle represents a sample, either bacterial or VLP. The bacterial and VLP communities from the same sample are connected by a line, where the red half of the line touches the VLP community, and the black half touches the bacterial community. The percent of total variation accounted for by each axis is shown in the axis label. The alignment of bacterial and VLP communities was highly significant (P = 0.004).

In order to test the effect of diet on the composition of the gut virome, we compared the distance between VLP communities in individuals both before and after they started their controlled diet (Fig. 5C). The distance between the gut viromes of individuals on the same diet (Fig. 5C, left) was significantly smaller at the end of their dietary treatment than it was at the start (P = 0.05, significance assessed by permutation of diet labels). There was no increase over time in virome similarity for individuals on different diets (Fig. 5C, right). A group of 39 VLP contigs were identified that changed in association with the dietary intervention (Supplemental Table S3). In comparison to the full data set of 7175 contigs, this subset showed a trend toward enrichment in Siphoviridae (P < 0.06) and depletion of Myoviridae (P < 0.08), which is suggestive of a phylogenetically distinct diet responsive bacteriophage population.

This dietary effect was also seen in a subject-by-subject analysis. For each subject, the distances between samples taken before the defined diet were compared, and distances between samples taken after defined diet were compared (Fig. 5D). As a contrast, distances between samples, where one was before-diet and the other was after-diet, were also compared. The distances were greatest for comparisons between before-diet and after-diet samples (Fig. 5D; P < 0.05, significance assessed by permutation of day labels). Thus, over the period of the dietary intervention, the viral community changed detectably to a new state.

The gut virome is dominated by prokaryotic viruses that prey on gut bacteria, raising the question of how changes in bacterial and viral communities are linked. We tested covariation by comparing 16S rDNA amplicon data from the gut DNA samples, which measured bacterial diversity with the VLP shotgun metagenomic data, which measured viral diversity (Fig. 5E). Variation was significantly correlated between bacterial and VLP communities (Mantel test, R = 0.40, P < 0.001). This correlation was significant when tested by subject label permutation (P < 0.001), indicating that interpersonal variation is correlated between these communities.

DNA isolation, PCR amplification, and purification

VLP DNA was isolated in the following manner (Fig. 1). Approximately 0.5 g of stool was homogenized in 40 mL of SM buffer (Sambrook and Russell 2001). Particulate matter was spun down at 4700g for 30 min and supernatant was filtered at 0.22 μm (PES filter, Nalgene). The filtrate was centrifuged on a CsCl step gradient (described in detail in Sambrook and Russell 2001; Thurber et al. 2009), and the 1.35–1.5 g/mL fraction was extracted from the column. We note that some viruses are known to have densities outside of this range (e.g., Tectiviridae, Poxviridae, Herpesviridae), and so would be lost during purification (Thurber et al. 2009); no attempt was made to isolate RNA viruses. Samples were treated with chloroform for 10 min, treated with DNase (Invitrogen) for 10 min at 37°C, and then extracted with the DNeasy Blood and Tissue Kit (Qiagen). VLP DNA was amplified with Genomiphi V2 polymerase (GE Healthcare).

Total stool DNA was isolated using the QIAamp DNA Stool Mini Kit (Qiagen), as described previously (Wu et al. 2010).

454/Roche sequencing methods

Both total DNA and amplified VLP DNA were randomly sheared and ligated to adapters using the 454 GS Titanium Rapid Library Preparation Kit with MID adaptors (454). Bacterial 16S amplicons were amplified with primers BSR 357-A (bar coded) and BSF8-B (not bar coded), which anneal to the V1–V2 region of the bacterial 16S rRNA gene, and include 454 sequencing adaptors (Wu et al. 2010). Samples were pooled and sequenced using GS FLX Titanium chemistry on a 454 Genome Sequencer (Margulies et al. 2005).

16S sequence analysis

Bacterial 16S sequences were analyzed using Pyronoise, assigned by comparison to the RDP database (Wang et al. 2007), and communities compared with each other using UniFrac (Lozupone and Knight 2005) scores, all within the Qiime software package (Caporaso et al. 2010). Sequences were filtered according to size (between 200 and 800 nt), ambiguous bases were removed (max = 2), homopolymer runs were truncated (max = 20), and primer mismatches removed (max = 0).

VLP contig sequence analysis

Functional comparison of VLP DNA and total DNA

Biochemical functionality was predicted according to similarity of ORFs to proteins within the COG database (Tatusov et al. 2003). Membership in a given COG category was determined by the proportion of reads (out of a given data set) that fell within ORFs that had significant similarity (E-value < 10⁻⁵) to a protein in the COG database in that category. Differences between VLP and total DNA were calculated via a two-sample t-test.

PHACCS analysis

Community structure was estimated by using PHACCS (v1.1.2) (Angly et al. 2005), implemented with the Octave (v3.2) environment using contig spectra as input, which were generated by Circonspect (v0.2.4). The contig spectra were generated by randomly sampling 10,000 sequences truncated to 100 bp in length, assembling at 98% identity and 35 bp overlap, and counting the number of contigs with one member, two members, etc. PHACCS was run on those spectra using a genome size of 50 kb, under a power law scenario.

We thank members of the Wu, Lewis, and Bushman laboratories for help and suggestions. This work was supported by the Human Microbiome Roadmap Demonstration Project UH2DK083981 (G.D.W., F.D.B., J.D.L., Co-PIs). We also acknowledge the Penn Genome Frontiers Institute and a grant from the Pennsylvania Department of Health (the Department of Health specifically disclaims responsibility for any analyses, interpretations, or conclusions); NIH AI39368 (G.D.W.); the Molecular Biology Core of The Center for Molecular Studies in Digestive and Liver Diseases (P30 DK050306); and The Joint Penn-CHOP Center for Digestive, Liver, and Pancreatic Medicine. We also acknowledge NIH instrument grant S10RR024525 and NIH CTSA grant UL1RR024134 from the National Center for Research Resources, and the Crohn's and Colitis Foundation of America. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Center for Research Resources or the National Institutes of Health.