Taxonomic and Functional Dysregulation in Salivary Microbiomes During Oral Carcinogenesis

Abstract

Exploring microbial community compositions in humans with healthy versus diseased states is crucial to understand the microbe-host interplay associated with the disease progression. Although the relationship between oral cancer and microbiome was previously established, it remained controversial, and yet the ecological characteristics and their responses to oral carcinogenesis have not been well studied. Here, using the bacterial 16S rRNA gene amplicon sequencing along with the in silico function analysis by PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States 2), we systematically characterized the compositions and the ecological drivers of saliva microbiome in the cohorts of orally healthy, non-recurrent oral verrucous hyperplasia (a pre-cancer lesion), and oral verrucous hyperplasia-associated oral cancer at taxonomic and function levels, and compared them with the re-analysis of publicly available datasets. Diversity analyses showed that microbiome dysbiosis in saliva was significantly linked to oral health status. As oral health deteriorated, the number of core species declined, and metabolic pathways predicted by PICRUSt2 were dysregulated. Partitioned beta-diversity revealed an extremely high species turnover but low function turnover. Functional beta-diversity in saliva microbiome shifted from turnover to nestedness during oral carcinogenesis, which was not observed at taxonomic levels. Correspondingly, the quantitative analysis of stochasticity ratios showed that drivers of microbial composition and functional gene content of saliva microbiomes were primarily governed by the stochastic processes, yet the driver of functional gene content shifted toward deterministic processes as oral cancer developed. Re-analysis of publicly accessible datasets supported not only the distinctive family taxa of Veillonellaceae and Actinomycetaceae present in normal cohorts but also that Flavobacteriaceae and Peptostreptococcaceae as well as the dysregulated metabolic pathways of nucleotides, amino acids, fatty acids, and cell structure were related to oral cancer. Using predicted functional profiles to elucidate the correlations to the oral health status shows superior performance than using taxonomic data among different studies. These findings advance our understanding of the oral ecosystem in relation to oral carcinogenesis and provide a new direction to the development of microbiome-based tools to study the interplay of the oral microbiome, metabolites, and host health.

Keywords: machine learning; microbiome dysbiosis; oral cancer; oral verrucous hyperplasia; saliva.

Figure 1

Differences in oral microbiomes among normal, OVH, and OSCC cohorts. (A, B) Adonis analysis based on (A) unweighted and (B) weighted UniFrac distance metrics shows the effect (R² ) of factors with the oral microbiome. * indicates FDR-adjusted p < 0.05 and ** indicates FDR-adjusted p < 0.01. (C, D) Principal coordinate analysis (PCoA) plots of taxonomic profiles based on (C) unweighted and (D) weighted UniFrac distance metrics. Marginal kernel densities visualize the distribution of microbial diversity along both axes. The pairwise PERMDISP reveals the dispersion effect (FDR-adjusted p < 0.05) between normal and OVH cohorts.

Figure 2

Core microbiome analysis. (A) Venn diagram of core microbiomes among cohorts. The core is defined as the species taxa present in saliva with ≥ 75% prevalence. (B) The fraction of core species number to overall species richness in each cohort. (C) LEfSe reveals the distribution of core species displaying the abundance significantly higher (LDA > log10³) among cohorts. The asterisk (*) indicates a taxon that was annotated only to the genus level. (D) Same as (C) but at the family level.

Figure 3

Distribution of signature pathways. The signature pathways, which abundances are significantly higher concerning each studied cohort, are detected using LEfSe. The inferred pathways are collapsed to each category based on Metacyc’s pathway ontology. Colored boxes indicate a higher rank of the categories.

Figure 4

Multiple-site beta diversity (Sørensen dissimilarity) and corresponding nestedness and turnover components. The dissimilarities were analyzed in terms of species and metabolic pathway profiles in each cohort.

Figure 5

Boxplots illustrate the null-model-based stochastic ratio of microbial taxonomic composition and functional profile based on Bray–Curtis dissimilarities. The simulated procedure was repeated 999 times with proportional occurrence frequency and richness.

Figure 6

Comparison of diversity and core analysis between taxonomic and functional profiles from this study and previous studies (Wolf et al., 2017; Zhao et al., 2017). (A–D) Venn diagrams reveal common core species/pathways (prevalence > 75%) in normal and OSCC cohorts, respectively. (Taxonomic profiles of normal (A) and OSCC (B) cohorts; metabolic pathway profiles of normal (C) and OSCC (D) cohorts).

Figure 7

Evaluating functional profile as an alternative signature for OSCC detection using machine learning with the datasets of this study and previous studies (Wolf et al., 2017; Zhao et al., 2017). (A) The mean accuracy ratios of 100 iterations of the randomly split dataset (80% training and 20% testing) were based on taxonomic and functional profiles. The accuracy ratio is defined as the predicted accuracy to the accuracy of a random guess. Vertical bars indicate 95% confidence intervals. (B, C) The 2D-density plots of ROC curves from 100 iterations demonstrate a higher mean AUROC using (B) functional profiles to distinguish OSCC from normal cohorts compared to that using (C) taxonomic profiles.