Evidence of a robust resident bacteriophage population revealed through analysis of the human salivary virome

Abstract

Viruses are the most abundant known infectious agents on the planet and are significant drivers of diversity in a variety of ecosystems. Although there have been numerous studies of viral communities, few have focused on viruses within the indigenous human microbiota. We analyzed 2 267 695 virome reads from viral particles and compared them with 263 516 bacterial 16S rRNA gene sequences from the saliva of five healthy human subjects over a 2- to 3-month period, in order to improve our understanding of the role viruses have in the complex oral ecosystem. Our data reveal viral communities in human saliva dominated by bacteriophages whose constituents are temporally distinct. The preponderance of shared homologs between the salivary viral communities in two unrelated subjects in the same household suggests that environmental factors are determinants of community membership. When comparing salivary viromes to those from human stool and the respiratory tract, each group was distinct, further indicating that habitat is of substantial importance in shaping human viromes. Compared with coexisting bacteria, there was concordance among certain predicted host-virus pairings such as Veillonella and Streptococcus, whereas there was discordance among others such as Actinomyces. We identified 122 728 virulence factor homologs, suggesting that salivary viruses may serve as reservoirs for pathogenic gene function in the oral environment. That the vast majority of human oral viruses are bacteriophages whose putative gene function signifies some have a prominent role in lysogeny, suggests these viruses may have an important role in helping shape the microbial diversity in the human oral cavity.

Figure 1

Plots of shared homologs for intra-subject and inter-subject comparisons of viromes. Databases were created for the reads from each virome, and homologs between viromes were determined based on significant blastN hits (E-score <10⁻⁵). The number of significant hits per E-score is shown on the x-axis and E-scores are shown on the y-axis. For E-score values ⩾180 (the equivalent of 0), the proportion of significant hits is shown above the dashed line. (a–e) Intra-subject comparisons for subjects 1, 2, 3, 4, and 5, respectively. (f) Inter-subject comparisons. The insets with Venn diagrams show the overall percentage of shared homologs amongst each virome.

Figure 2

Analysis of metabolic potential and virulence factor homologs in viromes from each subject. (a) Percentage of viromes from each subject devoted to various putative metabolic categories based on reads with significant blastX homology to known entries in the SEED database (E-score <10⁻⁵). (b) Heatmap of virulence factor homologs present in viral contigs for each subject at all time points. Virulence factors were defined as those gene sequences that contribute substantially and may be either directly or indirectly involved in pathogenesis, and homologs were determined based on significant (E-score <10⁻⁵) blastX homology to the Virulence Factor Database (Yang et al., 2008). Putative functional categories are listed on the left.

Figure 3

Percentage of integrases in viral contigs and putative prophage assemblies from Subject no. 1. (a) Contigs with integrase homologs are shown in blue and contigs without integrase homologs are shown in red for each subject across all time points. Putative Veillonella phage assemblies from Subject no. 1 on day 1 (b) and day 30 (c). Putative virulence factor homologs virE and yadA are shown in blue, and all other open reading frames are shown in green. (b) Viral contig (27 429 nucleotides, 1421 reads, average coverage 22 × ) on day 1 and (c) viral contig (27 904 nucleotides, 1302 reads, average coverage 20 × ) on day 30. In panels b and c, the genome of V. dispar ATCC17748 is shown, and the proportion of virome reads mapping to different portions of the genome are shown in purple.

Figure 4

Taxonomic assignments and residual plots comparing viruses and their bacterial hosts for all subjects at all time points. Phylum-level taxonomic assignments for putative viral hosts based on blastX best hits of contigs against the NCBI NR database are shown in panel a, and assignments for bacteria based on 16S rRNA sequences are shown in panel b. Genus-level residual plots for taxonomic assignments comparing bacterial taxonomy with putative viral host taxonomy are shown in panel c. The dashed lines represent significant residuals with P-values <0.01.

Figure 5

Heatmap of taxonomic assignments based on blastX best hits for viral contigs (a) and principal-coordinates analysis of viruses based on blastX best hits for viral contigs (b). At each heatmap time point, values are normalized by the total number of viral contigs. Principal-coordinates analysis was performed on Bray–Curtis values for viruses at all time points for each subject. Blue represents Subject no. 1, orange represents Subject no. 2, magenta represents Subject no. 3, gray represents Subject no. 4 and black represents Subject no. 5. A pooled respiratory virome is represented by green and 11 individual stool viromes are represented by purple.