Core microbiome across dietary patterns

A vegan and omnivore core saliva microbiome was defined as the genera present in >95% of omnivores and vegans, respectively. Twenty-three genera present in both core microbiomes constituted a common saliva microbiome across dietary patterns accounting for 97.0 ± 2.2% (mean ± standard deviation) of all reads in each individual (Supplementary Fig. S1A). Compositionally, the core microbiota was dominated by members of the three major phyla Bacteroidetes, Firmicutes and Proteobacteria (Supplementary Fig. S1B), with Prevotella, Veillonella, Neisseria and Streptococcus as the predominant genera, each with average relative abundance >5%, albeit with substantial inter-individual variation.

We identified 12 operational taxonomic units (OTU) that were present in all individuals (Supplementary Table S2). Among these were OTUs assigned to Neisseria subflava, Haemophilus parainfluenzae, Prevotella melaninogenica, Veillonella dispar and Veillonella parvula, as well as unclassified Streptococcus spp., Granulicatella spp. and Campylobacter spp. Interestingly, these 12 core OTUs accounted for more than half of all reads (51.5 ± 7.7%).

Analyses of co-occurrence and co-exclusion revealed interesting patterns among the 23 core genera (Supplementary Fig. S2). A high degree of pairwise co-occurrence was observed between Porphyromonas spp. and Fusobacterium spp. (ρ = 0.74; Q < 10−15) and between Porphyromonas spp. and Neisseria spp. (ρ = 0.71; Q < 10−15). Conversely, pronounced co-exclusion was observed between Veillonella spp. and Porphyromonas spp. (ρ = −0.73; Q < 10−15), Veillonella spp. and Neisseria spp. (ρ = −0.75; Q < 10−15), and Neisseria spp. and Prevotella spp. (ρ = −0.70; Q < 10−15).

Alpha and beta diversity in vegans and omnivores

In the cohort as a whole, mean microbial richness in terms of observed and estimated (Chao1) number of OTUs was 201 ± 35.7 and 261 ± 51.0, respectively. We did not observe any difference in richness nor overall diversity as assessed by Shannon’s index and Simpson’s reciprocal index when contrasting vegans and omnivores (Supplementary Fig. S3). However, principal coordinate analysis (PCoA) revealed a subtle difference in beta-diversity in vegans and omnivores (Supplementary Fig. S4). Using permutational multivariate analysis of variance (PERMANOVA) to contrast dietary patterns, we found significant differences in community structure as assessed by Bray-Curtis (R2 = 2.1%, P = 0.008) and weighted UniFrac (R2 = 2.6%, P = 0.019) distances, but no difference for unweighted UniFrac (R2 = 0.8%, P = 0.125), indicating that the microbial communities in vegans and omnivores are phylogenetically similar, with difference in community structure driven by varying abundances of OTUs present in both vegans and omnivores.

Compositional differences in vegans and omnivores

For analysis of compositional differences between vegans and omnivores a subset of OTUs with a mean relative abundance >0.01% and prevalence >50% across diet groups were considered. Differential abundance was observed at all taxonomic levels below phylum level (Supplementary Fig. S5). A total of 22 OTUs were differentially abundant at an FDR < 10% and 10 OTUs were differentially abundant at an FDR < 5% (Fig. 1; Supplementary Table S3). Among the differentially abundant OTUs were oral and upper respiratory tract commensals like Neisseria subflava, Haemophilus parainfluenzae, Rothia mucilaginosa, and Capnocytophaga spp. which were more abundant in vegans, and Prevotella melaninogenica and Streptococcus spp. which were more abundant in omnivores. Campylobacter rectus and Porphyromonas endodontalis, species associated with periodontal disease, were also more abundant in vegans.

Figure 1 Differential abundance of operational taxonomic units. Volcano plot of estimated log 2 fold difference in operational taxonomic unit (OTU) abundance between vegans and omnivores and corresponding Benjamini-Hochberg adjusted P-values (Q) from negative binomial Wald tests as implemented in the DESeq2 R package. The red dotted line indicates the 10% false discovery threshold. Prevalence indicates percentage of participants in which a given OTU is present. Abundance indicates mean relative abundance (‰) of a given OTU. Name of OTUs differentially abundant at an FDR ≤ 5% are given at the lowest classified rank in Greengenes [Greengenes ID]. See Supplementary Table S2 for a full list of OTUs differentially abundant at an FDR < 10%. p, phylum, o, order. f, family. g, genus. s, species. Full size image

Salivatypes

Clustering of samples based on genus abundance, using a partitioning around medoids (PAM) approach as previously described13, identified two semi-separate clusters; cluster I characterized by higher abundance of Neisseria and Porphyromonas, and cluster II characterized by higher abundance of Prevotella and Veillonella (Fig. 2). Both clusters included vegan and omnivore samples, with 46.2% of vegans and 34.1% of omnivores assigned to cluster I, and 53.8% of vegans and 65.9% of omnivores assigned to cluster II (P = 0.15, Fisher’s exact test). To address the issue of salivatypes as biological gradients rather than discrete features, we plotted the ratio of the predominant genera, Neisseria and Prevotella, against the first principal coordinate axis of the Jensen-Shannon divergence used to define the salivatype clusters, showing the discrete salivatypes to be extremes of a continuum (Fig. 2B). Interestingly, the ratio of Neisseria to Prevotella was significantly higher in vegans (P = 0.0008, Wilcoxon rank-sum test), reflecting an association between dietary patterns and the compositional features underlying genus based clusters.

Figure 2 Salivatypes in vegans and omnivores. (A) Principal coordinates analysis visualizing salivatype clusters based on partitioning around medoids of Jensen-Shannon distance (JSD). Ellipses cover 67% samples in each cluster. (B) The ratio of Neisseria to Prevotella along the first principal coordinate axis of Jensen-Shannon distances used to build the salivatype clusters. (C–F) Relative abundance of the main contributors to each salivatype cluster. Differential abundance of each genus on which the clusters were build was tested using a Wilcoxon rank-sum test and the genera with P values < 10−10 are depicted. Boxes represent interquartile range (IQR), with the inside line representing the median. Whiskers represent values within 1.5 × IQR of the first and third quartiles. Circles represent outliers beyond the whiskers. Full size image

Association of nutrients with microbial composition in saliva

In order to assess the impact of collinear nutrients on the diversity of the oral microbial community, we tested the association between four measures of alpha-diversity and the first 10 principal components (PC) of the daily macro- and micronutrient intake, capturing >75% of the variation in daily nutrient intake. Using multiple regression including all 10 PCs simultaneously, we found that PC3 of the diet data was negatively associated with both Shannon’s index (0.7% decrease per unit increase in PC3; 95%CI: 1.2–0.2%; Q = 0.05) and Simpson’s reciprocal index (2.6% decrease per unit increase in PC3; 4.4–0.8%; Q = 0.09). There was no association between PC3 and neither observed nor estimated OTU count, indicating an effect of the principal component on the evenness of the community rather than the richness.

Using coinertia analysis to assess the overall congruence between variation in the salivary microbiota composition at OTU level and the variation in average daily intake of macro- and micronutrients (Supplementary Fig. S6), we found a moderate but significant correlation (RV = 0.08, P = 0.002).

We also tested the impact of collinear nutrients on microbial community structure and community membership using PERMANOVA, showing multivariate associations (the effect of a given PC adjusted for the remaining 9 PCs) between three diet PCs and three measures of beta diversity (Fig. 3). PC3 explained 5.2% (Q = 0.05) and 2.1% (Q = 0.06) of the variation in weighted UniFrac distance and Bray-Curtis dissimilarity, respectively. Similarly, PC2 explained 3.0% (Q = 0.05) and 2.2% (Q = 0.05) of the variation in weighted UniFrac and Bray-Curtis distances, respectively. PC5 was associated with weighted UniFrac (R2 = 3.6%, Q = 0.06), but not with Bray-Curtis dissimilarity. PC2 was negatively associated with dietary fibre (contributing >4% to the PC; Supplementary Fig. S7A). PC3 was negatively associated with the medium-chain fatty acids (MCFA) octanoic (caprylic, C8:0), decanoic (capric, C10:0), and dodecanoic (lauric, C12:0) acid (each contributing >7.5% to the PC; Supplementary Fig. S7B), all three constituents of coconut oil and palm kernel oil. PC5 was negatively associated with starch, and positively associated with the omega-3 polyunsaturated fatty acids (PUFA) stearidonic acid (SDA, C18:4n3), eicosapentaenoic acid (EPA, C20:5n3), docosapentaenoic acid (DPA, C22:5n3) and docosahexaenoic acid (DHA, C22:6n3) as well as the omega-9 and omega-11 mono-unsaturated fatty acids (MUFA) nervonic acid (C24:1n9) and cetoleic acid (C22:1n11), all predominantly reflecting intake of fish (Supplementary Fig. S7C). Combined, the 10 first dietary PCs explained 18.7%, 12.5%, and 8.1% of the variation in weighted UniFrac distance, Bray-Curtis dissimilarity, and unweighted UniFrac distances, respectively.

Figure 3 Effect of diet principal components on alpha and beta diversity. (A) Forest plot of effect sizes and corresponding 95% confidence intervals of the first ten diet principal components (PC) on observed richness, estimated (Chao1) richness, Shannon’s diversity index, and Simpson’s reciprocal index. Associations were tested using multiple regression including all ten PCs simultaneously. (B) Heatmap of variance in beta diversity explained by each of the ten first diet PCs as estimated by permutational analysis of variance. Q values are given for associations significant at an FDR ≤ 10% only. Full size image

In order to clarify the compositional changes underlying the association between diet and alpha and beta diversity we focused on the subset of OTUs with a mean relative abundance >0.01% and prevalence >50% across diet groups, divided samples into quartiles according to each of PC2, PC3 and PC5, and contrasted the first and fourth quartile for each of the principal components. For each of the three PCs we found several associations significant at an FDR < 5% with substantial effect sizes ranging from two-fold to four-fold differences (Supplementary Table S4). PC2 (low fibre intake) was negatively associated with two OTUs assigned to genus Capnocytophaga and two OTUs assigned to genus Neisseria, one of which was annotated at species level as Neisseria subflava. Another OTU assigned to Neisseria subflava was negatively associated with PC5 (high piscine PUFAs and low starch). PC3 (low caprylic, capric, and lauric acid intake) was negatively associated with two OTUs assigned to genus Prevotella, one OTU assigned to genus Leptotrichia and one OTU assigned to genus Selenomonas. One OTU assigned to species Veillonella dispar was positively associated with PC2, whereas another OTU also assigned to V. dispar was negatively associated with PC5. Similarly, one OTU assigned to family Neisseriaceae was positively associated with PC3, while another OTU assigned to Neisseriaceae was negatively associated with PC5. PC5 was also positively associated with one OTU assigned to genus Actinobacillus, one assigned to genus Lautropia, and one assigned to genus Porhyromonas.

Association of microbial composition with inflammatory biomarkers

In order to identify salivary microbiota associated with low-grade inflammation we divided samples into quartiles according to serum CRP concentration and white blood cell (WBC) count and contrasted the upper and lower quartile for each of the principal components, focusing again on the subset of OTUs with a mean relative abundance >0.01% and prevalence >50% across diet groups. We identified six OTUs associated with inflammatory makers (Supplementary Table S5); five were associated with CRP and one with WBC. One OTU assigned to Haemophilus parainfluenzae was associated with both CRP and WBC, but in opposite direction (positively associated with WBC and negatively associated with CRP).

Functional differences in vegans and omnivores

Based on rarefied OTU counts the functional potential of the microbial communities in vegans and omnivores was predicted using PICRUSt14. A mean weighted Nearest Sequenced Taxon Index (weighted NSTI) score of 0.03 ± 0.02 indicated high reference genome coverage. Considering the resulting 6,909 KEGG orthologous groups (KOs), we did not find a difference in functional richness between vegans and omnivores. However, at an FDR < 0.10 we did find significantly higher alpha diversity indices (Shannon’s index and Simpson’s reciprocal index) in vegans compared to omnivores, indicating that the difference in functional alpha diversity is driven by difference in abundance of shared features rather than presence or absence of specific KOs. PCoA ordination of Bray-Curtis dissimilarity indicated a subtle (R2 = 4.0%) but significant (P = 0.002) difference in functional beta diversity (Supplementary Fig. S8). In fact, diet explained a larger proportion of the variation in the genomic potential than it did variation in community structure (4.0% vs. 2.1% variation in Bray-Curtis dissimilarity explained, respectively).

When collapsing KOs into KEGG pathways we identified 183 pathways within the overall topics of metabolism, environmental information processing, genetic information processing and cellular processes that were present in more than 10% of all individuals, 85 of which were differentially abundant between vegans and omnivores at an FDR < 5%; 46 pathways were enriched in vegans and 39 enriched in omnivores (Fig. 4). In general, effect estimates were small. Of the pathways enriched in vegans, 50% had a log 2 fold-difference below 0.06, corresponding to less than a 4% increase compared to omnivores. Of the pathways enriched in omnivores, 50% had a log 2 fold-difference above −0.04, corresponding to less than a 2.5% increase compared to vegans. Additionally, a substantial proportion of the differentially abundant pathways constituted a very low proportion of the genomic potential within each sample; 23 pathways (27%) covered less than 1‰ each.