Signatures of Mucosal Microbiome in Oral Squamous Cell Carcinoma Identified Using a Random Forest Model

Introduction

Oral cancer is a fatal malignancy and the sixth most common cancer worldwide. Oral squamous cell carcinoma (OSCC) is the most frequently occurring oral cancer, accounting for over 90% of all oral cancers.¹ In recent years, the incidence of OSCC has shown an increasing trend, which constitutes a major global health problem. OSCC is associated with a high rate of relapse, unpredictable metastasis and oral maxillofacial damage.^2,3 The 5-year survival rate of patients with OSCC is approximately 50–60%.³

Mounting evidence supports the association between oral microbiome and OSCC. It has been suggested that the oral microbiome plays a role in the development of oral cancer.^4–6 There are three potential mechanisms by which the oral microbiome is involved in the pathogenesis.⁶ Firstly, oral bacteria can induce chronic inflammation. Inflammatory mediators produced from the oral mucosa promote cell hyper-proliferation, mutagenesis, oncogene activation and angiogenesis, leading to carcinogenesis. Secondly, certain oral bacteria are directly involved in the development of OSCC through enhancing cell proliferation, causing cytoskeletal rearrangement and activation of NF-κB pathways, and inhibiting cellular apoptosis. Finally, some metabolites from oral bacteria were shown to possess carcinogenic potentials.⁶

A number of bacteria from the oral microbiome have been shown to be associated with OSCC. The relative abundance of Fusobacterium sp., Porphyromonas gingivalis, Streptococcus sp., Peptostreptococcus sp., and Prevotella sp. are increased in OSCC.^7,8 Zhao et al reported significant enrichment of Fusobacterium, Dialister, Peptostreptococcus, Filifactor, Peptococcus, Catonella, and Parvimonas in OSCC patients compared with matched controls using linear discriminant analysis (LDA) effect size (LEfSe) analysis.⁹ In a prospective study, involving examination of mouth swab microbiome, Parvimonas micra and Neisseria sicca were associated with reduced risk of oral cavity cancer.¹⁰ The identified OSCC associated oral bacteria, however, vary among studies. This is partially attributable to the variations between studies with respect to the sample types used for the analysis of oral microbiome. There are substantial variations in the composition of the oral microbiome between various sample types including tissue, saliva, and mouth swab.¹¹ Few studies have characterized the OSCC related oral microbiome from tissue, which is closely associated with oral mucosa.

Features of the oral microbiome associated with OSCC have been explored by comparing OSCC patients with healthy controls, or by comparing tumor sites with the surrounding normal tissue.^9,10 It has been alerted microbial diversity and a distinctive profile of the oral microbiome in OSCC.^8–10 Many bacteria of the oral microbiome are enriched or depleted in OSCC.^8,9 Recent studies have also revealed that the collective behaviors of the microbial community are more closely related to the development of cancer in comparison with individual bacteria genera.¹² To further explore how the oral microbiome is involved in the carcinogenic process, it is of importance to identifying signatures discriminative for the OSCC associated oral microbiome. The random forest (RF) method is a machine learning program capable of identifying an optimal set of variables with high discriminative power from a large number of dependent or independent variables. It is known for its ability to render complex dependency patterns between the outcome and the covariates.¹³ The method has been widely used to identify microbial signature for colorectal carcinoma, gastric cancer, and pancreatic cancer.^14–16 In this study, we applied the RF model to explore the signatures of the OSCC associated oral microbiome and validated these in an independent cohort.

Patients and Methods

Patients and Sample Collection

A total of 24 patients at different stages of OSCC were enrolled in the study (Table 1). The OSCC stages were classified according to the American Joint Committee on Cancer (AJCC) staging manual (8th edition).¹⁷ Of these patients, 19 were male. The average age of patients was 61.1 ± 12.4 years. None of the enrolled subjects had history of betel chewing, or history of diabetes mellitus or other severe complications including heart, liver, or renal failures. These patients had not taken antibiotics for at least two weeks immediately prior to sampling. Tissues from the site of tumor and paracancerous tissue from 2 cm outside the edge of the tumor were sampled. The paracancerous tissues are histologically normal. All samples were stored at –80°C until further processing. This study was approved by the Research Ethics Committee of Qingdao Municipal Hospital, China. Prior to sample collection, written informed consent for the study was obtained from all participants. All experiments were performed in accordance with approved guidelines and regulations of the institution.

Table 1 Clinical Characteristics of 24 Patients with Oral Squamous Cell Carcinoma

Analyses of the Oral Microbiome

To analyze the microbial communities of the oral mucosa, genomic DNA was extracted from the oral mucosa samples as previously reported.¹⁸ The variable V3–V4 region of the 16S rRNA gene was PCR amplified with primers 341F/806R to generate the amplicon libraries. Sequencing was performed on an Illumina PE250 platform (Illumina, Hayward, CA, USA). A total of 16,333,787 paired-ends reads were obtained. After quality control and filtration, 2,163,267 reads were produced with an average of 45,068 reads per sample. The reads were analyzed using UPARSE.¹⁹ Following global trimming at 250 nucleotides, reads were dereplicated, and singletons were discarded. Subsequently, reads were clustered into operational taxonomic units (OTUs) assuming 97% identity. Chimeric reads were then removed. Taxonomy assignment was performed using UClust.²⁰ Analyses of alpha and beta microbial diversity were conducted as described previously.¹⁸ Comparisons of the relative abundance of taxa between groups were performed using version 1.0 of LEfSe.²¹ A LDA value greater than 3 with a p value less than 0.05 was considered to be statistically significant. To analyze the correlation network, Spearman correlation coefficients were computed between the genera. Correlations that had an absolute Spearman coefficient value greater than or equal to 0.6 with a p value less than 0.05 were transformed into links between two genera in the genus network. Cytoscape v3.7.1 was then used to construct network figures. PICRUSt (v1.1.1) was used to predict the functions of the microbial community according to the reference data set of GreenGenes 16S rRNA database (v13.5).²² The accuracy of the predicted functions was assessed by the nearest sequenced taxon index (NSTI). The predicted functions were categorized with KEGG orthology. STAMP (v2.1.3) was used to compare functional differences between cancerous and paracancerous groups.²³

Statistical Analyses

To identify microbial signatures capable of distinguishing cancerous lesions from paracancerous lesions, a RF model was built using the AUC-RF algorithm.²⁴ The input variables comprised the relative abundance of taxa. A taxon was included only if it was present in more than 20% samples and had a relative abundance of greater than 0.05%. A 20-times repeated 10-fold cross-validation of the RF model was performed. The performance of the RF model was demonstrated by receiver operating characteristic (ROC) curve analysis.²⁵ Validation of OSCC associated bacteria was explored in an independent cohort in Shanghai (China) with 80 samples, of these, 40 samples were from OSCC surface lesion and 40 samples were from anatomically matched normal sites (BioProject accession number: PRJNA362794).⁹

The Mann–Whitney U-test was performed to detect significant differences in alpha diversity or relative abundance between disease groups. Significant P values were corrected for multiple testing using the Benjamini‐Hochberg procedure, with a cut-off false discovery rate (FDR) – adjusted P value lower than 0.1.²⁶ Constrained PCoA is a flexible method which allows for constrained ordination on the basis of any distance or dissimilarity measure.²⁷ It was performed in R to explore the influence of clinical stage on the microbiome of OSCC.

Data Availability

The next-generation sequencing datasets have been submitted to NCBI (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA597251/). The BioProject accession number is PRJNA597251.

Results

Differences in Microbial Community Between Cancerous and Paracancerous Tissues

The microbial community between the cancerous lesions (C) and matched paracancerous tissues (P) was compared in 24 patients. Both Chao1 and Shannon indices which estimate the biodiversity of microbiome, were not significantly different between group C and group P (P > 0.05 for both) (Figure 1A and B). On multiple linear regression analyses, the alpha diversity indices of the microbiome showed no association with age, sex, smoking, alcohol consumption, degree of differentiation, or clinical stage in either group C or group P. Principal coordination analyses (PCoA) showed that the microbial community structure of group C was distinct from that of group P (P = 0.048) (Figure 1C). Constrained PCoA analyses were performed to explore the influence of clinical stage on the OSCC microbiome. The results demonstrated an apparent separation between the clinical stage I, II, and III in group C, especially between stages II and III. This suggested the progression of OSCC was associated with alteration of oral mucosa microbiome (Figure 1D). In group P, however, the microbial community structure was not discernible in the plots between stages (Figure 1E).

Figure 1 Profiles of the oral microbiome in OSCC. The alpha diversity was estimated using Chao1 (A) and Shannon (B) indices. PCoA plots for comparing community structure between paracancerous and cancerous tissues (C), between clinical stages in cancerous tissues (D) and paracancerous tissues (E). The composition of microbiome at the phylum level was compared between paracancerous and cancerous tissues (F). (P) paracancerous tissues; (C) cancerous tissues.

Altered Composition of the Microbiome in OSCC

Compositional analyses showed the phyla of Proteobacteria, Firmicutes, Bacteroidetes and Fusobacteria were dominant in the microbiome of both groups (Figure 1F). The relative abundance of Firmicutes, Bacteroidetes, Fusobacteria, and Spirochaetes was increased in cancerous lesions. Differences with respect to the relative abundance of Fusobacteria (9.80% vs 1.70%) and Spirochaetes (3.63% vs 0.18%) were statistically significant using Mann–Whitney U-test (P = 0.043 and 0.032, FDR-adjusted P = 0.089 and 0.083, respectively). The abundance of Proteobacteria, Actinobacteria, and Cyanobacteria in cancerous lesions was lower compared to paracancerous tissues. The differences of relative abundance of Actinobacteria (0.77% vs 1.78%) and Cyanobacteria (0.23% vs 1.33%) were statistically significant (P = 0.021 and P = 0.013, FDR-adjusted P = 0.067 and 0.051, respectively). At the genus level, LEfSe analyses identified 15 genera which had a LDA score >3.0 (Figure 2A). The relative abundance of Fusobacterium, Treponema, Streptococcus, Peptostreptococcus, Carnobacterium, Tannerella, Parvimonas, and Filifactor were increased in cancerous lesions (Figure 2B). Of these, Fusobacterium had the highest LDA score. In contrast, the relative abundance of Streptophyta, Brevundimonas, Paenibacillus, Microbacterium, Desulfovibrio, Mucispirillum, and Arthrobacter were decreased in cancerous lesions (Figure 2C).

Figure 2 LEfSe analyses of microbiome composition between paracancerous and cancerous tissues. Bacteria genera enriched in cancerous tissues had a positive LDA score, while those depleted had a negative score. Bacteria with LDA scores >3 are shown (A). Relative abundance of genera enriched in cancerous tissues (B). Relative abundance of genera enriched in paracancerous tissues (C). *P < 0.05; **P < 0.01; ***P < 0.001. All FDR-adjusted P values were <0.05.

Microbial Signatures Associated with OSCC

To determine microbial signatures capable of discriminating cancerous lesions from paracancerous tissues, an RF classifying model was built with the AUC-RF algorithm. The results showed a minimal set of 12 bacteria genera that maximally differentiated cancer sites from paracancerous sites. It contained Microbacterium, Streptococcus, Brevundimonas, Peptostreptococcus, Streptophyta, Paenibacillus, Actinomyces, Parvimonas, Mucispirillum, Arthrobacter, Rhizobium and Fusobacterium (Figure 3A and B). Most of these (10/12) have been identified to have a LDA value > 3.0 by the LEfSe analyses, including four genera (Fusobacterium, Parvimonas, Peptostreptococcus and Streptococcus) enriched in cancerous lesions and six depleted genera (Arthrobacter, Brevundimonas, Microbacterium, Mucispirillum, Paenibacillus, and Streptophyta). RF models trained with this optimal set of features resulted in an out-of-bag error rate of 33.33%. To assess the model classifying accuracy, a 20-times repeated 10-fold cross-validation was performed. The AUC value from cross-validations was 0.82 [95% confidence interval (CI): 0.70–0.94] (Figure 3C). This set of 12 genera was further validated in an independent Chinese Shanghai cohort with AUC of 0.69 (95% CI: 0.54–0.80) (Figure 3D). The results suggested that these oral mucosal bacteria associated with OSCC may predict the risk of OSCC.

Figure 3 Identification of microbial signature associated with OSCC. The random forest model was constructed using AUC-RF algorithm based on bacteria which were present in more than 20% samples and had a relative abundance of over 0.05%. The model containing 12 genera was selected as the optimal model based on the highest OOB-AUC value obtained from the backward elimination process performed using the AUC-RF algorithm with the median decrease Gini (MDG) importance measure (A). MDG of selected genera in the optimal set (B). The ROC curve of the optimal model for classifying cancerous tissues from paracancerous tissues (C). The ROC curve of validation in a independent cohort in Shanghai (D).

Co-Occurrence Network in OSCC

The co-occurrence network differed substantially between cancerous lesions and paracancerous tissues (Figure 4). The network was simplified in group C with a decreased number of nodes (Table 2). The network density and average degree were also lower in group C, which demonstrated a reduction in the network complexity. The clustering coefficient of the network was lower, leading to an increased number of components and more fragmented sub-networks in group C (Figure 4).

Table 2 Features of Co-Occurrence Network of the Oral Microbiome

Figure 4 Correlation network in paracancerous tissues (A) and cancerous tissues (B). The correlation coefficient was calculated with Spearman’s rank correlation test (|r| ≥0.6). Cytoscape was used for network construction.

Predicted Functions of the OSCC Microbiome

Functional capacities of the OSCC microbiome were predicted using PICRUSt. Compared with paracancerous tissues, the average relative frequencies of pathways related to metabolism, especially amino acid metabolism, were decreased at KEGG levels 1 and 2 (Figure 5A and B). Furthermore, reduced relative frequencies in OSCC were found for pathways of KEGG level 3 related to amino acid metabolism including phenylalanine metabolism, tryptophan metabolism, arginine and proline metabolism, valine, leucine and isoleucine degradation, lysine degradation and tyrosine metabolism (Figure 5C). Nevertheless, the relative frequencies of methane metabolism, glucose-related metabolism, such as phosphotransferase system (PTS) and glycolysis, were significantly enriched in cancer lesions (Figure 5D). These alterations suggested functional disorders of amino acid metabolism and increased capacity of glucose utilization in the microbiome of OSCC. To validate these results, the predicted functions were compared between two randomized groups. The results showed no significant differences between two groups at KEGG levels 1–3 (Supplementary Table 1).

Figure 5 Differential functions predicted using PICRUSt between paracancerous and cancerous tissues. The function of metabolism at level 1 in KEGG pathway (A) and amino acid metabolism at level 2 (B). The decreased (C) or increased (D) relative frequency of functions of amino acid metabolism at level 3 in cancerous tissues. Differences between groups in the predicted functions were compared using STAMP. Statistical differences are considered when P < 0.05.

Discussion

Our results demonstrated significant differences in the β-diversity (including PCoA and constrained PCoA analysis) of the oral microbiome between cancerous and paracancerous tissues. These findings indicate altered community structure of the oral microbiome in cancerous lesions. However, the α-diversity including Chao1 and Shannon indices showed no significant differences. Consistent with our results, previous studies have found no significant differences in α-diversity between cancerous lesions and paratumoral mucosa in patients with colorectal cancer.²⁸ Constrained PCoA analyses demonstrated an apparent separation among the different clinical stages in cancerous lesions of OSCC. It should be noted that the sample size in different stages was small. Thus, the observed association between variations of the oral microbiome and cancer progression needs to be demonstrated in studies with a large sample size. Nonetheless, in a previous study, the oral microbiome in patients with stage IV OSCC was different from that in healthy controls.²⁹ The microbiome communities showed a progressive shift with the progression of OSCC stage.²⁹ These findings are consistent with our results wherein the community structure of the oral microbiome was found to changes with the stages of OSCC, indicating the potential involvement of the oral microbiome in cancer progression.

In recent studies, microbial signatures identified using RF machine learning algorithm were used to predict different diseases.^14–16 The identified signatures can be viewed as a minimized representation of the microbial community for exploring mechanisms of microbial dysbiosis and targets for manipulating the dysbiotic microbiome. In a study by Zackular et al, the final numbers of colon tumor were quantitatively predicted on the basis of the initial composition of the gut microbiome using RF algorithm.³⁰ In another report, a minimal set of 22 viral genera were identified that differentiated CRC patients from controls and predicted CRC outcomes.¹⁴ An RF model for predicting exacerbation of radiation-induced mucositis was generated using oral microbiome, it showed a high predictive accuracy (AUC = 0.89).³¹ Currently, few studies have been conducted to identify the microbial signature in OSCC. Our study revealed an RF model consisting of 12 relevant bacteria with differential abundance between the cancerous lesions and paracancerous tissues of OSCC. The trained model was capable of predicting OSCC with considerably high accuracy (AUC = 0.82). Moreover, the validation in an independent cohort showed a moderate accuracy. Monitoring the composition and abundance of these 12 bacteria may help predict the risk of OSCC, facilitate early diagnosis and treatment of OSCC and help improve the prognosis of OSCC.

In our study, the LDA scores of 15 genera were >3.0 when comparing the oral microbiome of cancerous lesions with that of paracancerous tissues. Furthermore, these bacteria showed differential abundance between the two groups, suggesting a close association between these bacteria and OSCC. The majority (10/12) of bacteria in the optimal set identified by RF analyses were contained in these 15 bacteria genera. The highly consistent results between these analyses demonstrate that variations of the optimal set of bacteria reflect the dysbiotic features of the oral microbiome associated with OSCC, accounting for remarkable capacity of the set to discriminating OSCC from normal tissues. Of the optimal set of bacteria, Fusobacterium, Parvimonas, Peptostreptococcus and Streptococcus were enriched in cancerous lesions. The association of these four bacteria with OSCC has been reported. It has been found an increase in the abundance of Fusobacteria in OSCC, especially in the late stage.^11,32 A significant enrichment of Fusobacterium, Peptostreptococcus, Peptococcus and Parvimonas in OSCC surface lesions has been demonstrated.^9,33 Variations in the abundance of Peptostreptococcus and Parvimonas in the salivary microbiome have been observed between normal individuals, oral precursor lesion and cancer patients.³⁴ Infection with S. anginosus occurs frequently in OSCC.³⁵ Species of these genera have been isolated from the oral tumor and metastatic lymphatic nodes.^36,37 These results support a close association between these enriched bacteria and OSCC. In addition, these enriched bacteria may have carcinogenic potential. Previous studies have shown that F. nucleatum may induce tumorigenesis through the inflammatory NF-κB signaling pathway.^38,39 F. nucleatum expresses the adhesion molecule FadA which binds to E-cadherin on epithelial cells, inhibiting its tumor-suppressor activity.⁴⁰ Streptococcus spp. possesses alcohol dehydrogenase activity leading to production of the carcinogenic acetaldehyde, which is potentially involved in the occurrence of oral cancer.⁴¹ These findings suggest that our RF analyses identified a group of bacteria with carcinogenic potentials enriched in OSCC. In addition to these bacteria, we also found depletion of six bacteria genera in OSCC. Further studies are required to unravel the causes of variations in the abundance of these bacteria and their roles in carcinogenesis.

Bacterial interaction is a determinant of the homeostasis of the microbiome. Studies have shown reduced network complexity in peritumoral or tumoral tissues of the stomach in comparison with normal gastric mucosa.⁴² In our study, strongly interacted genera were picked up according to |r| ≥0.6 to construct microbial network. Our results demonstrated a simplified microbial network in cancerous lesions of OSCC, with the reduction in the number of nodes, network density, average degree, and clustering coefficient in cancer sites, suggesting fragmented and simplified network. These findings indicated the loss of certain bacterial interactions in OSCC, leading to disruption of the homeostasis of oral microbiome. However, a recent study found increased network complexity in OSCC.⁹ The contradictory results may be attributed to the difference in the cut-off value of correlation coefficient used bacterial species for the construction of microbial network.

The functional profile in the cancer-associated microbiome shows substantial changes. In a previous study, the functions related to protein and amino acid metabolism (such as valine, leucine and isoleucine, phenylalanine, tyrosine and tryptophan biosynthesis) showed an inverse association with OSCC progression.²⁹ Similar results were obtained in a study of saliva metabolomics wherein phenylalanine metabolism was found down-regulated in OSCC.⁴³ Consistent with this, in our study, the functional prediction of oral bacterial communities showed a significant decrease in genes involved in amino acid metabolism (such as phenylalanine metabolism, tryptophan metabolism, arginine, and proline metabolism) in OSCC lesions. However, methane metabolism and glucose-related metabolism (such as PTS and glycolysis) were significantly enriched in cancer lesions. Tryptophan metabolism is associated with the evasion of immune surveillance by cancer cells.⁴⁴ Decreased levels of tryptophan were shown to promote the generation of immune-suppressive T cells promoting cancer development.⁴⁵ Cancer cells appear to utilize glycolytic metabolism for energy production and precursors for the increased rate of anabolic processes.⁴⁶ These suggest the gene functional disorders of amino acid metabolism and increased capacity of glucose utilization in the microbiome are associated with OSCC.

In summary, this study revealed dysbiosis in OSCC featured with functional and compositional changes. We identified a microbial signature comprising 12 bacteria associated with OSCC, which may be used as a biomarker for clinical prediction of the risk of OSCC. Further studies are required to explore the potential clinical application of the microbial signature.