Clinical characteristics of the cohortsFresh frozen primary tumors from 50 patients with ACC, stage 1–4, and no prior therapy at the time of surgery were subjected to 16S RNA gene sequencing to characterize the intratumoral bacterial composition, diversity, structure, and association with clinicopathological factors (Supplementary Fig. S1). Adjacent nontumoral tissue from thirty-three patients was also sequenced for comparative analysis of microbial communities with those in the tumor. The tumors of 44 patients (88%) were analyzed using RNA-seq in our previous study4. A classification of these tumors into 2 molecular subtypes, referred as ACC-I and ACC-II, was available for correlative analysis with the intratumoral microbiome. Patients with both molecular subtypes (45% ACC-I and 55% ACC-II) were included in the study cohort (Table 1, Supplementary Table). The primary sites of ACC in the cohort were glands characterized by the following locations: maxillary sinus, base tongue, palate, parotid, trachea, sublingual, lacrimal, and submandibular. The sites are close to the oral, throat, nasal, and eye ocular surface body sites characterized by the Human Microbiome Project. Most of the patients were males (66%), had perineural invasion (PNI) (84%), and stage 3–4 tumors at diagnosis (70%) with solid or cribriform histology (84%). All the patients underwent primary tumor resection, and most of them were treated with postoperative RT (88%) and/or chemotherapy (38%).Table 1 Clinicopathological characteristics of the ACC cohort.The taxonomic structure of the bacterial community in ACC tumors and adjacent normal tissue is associated with bacterial richness16S RNA sequencing of 50 ACC tumors revealed a median number of 200 putative bacterial species per sample, ranging from 52 to 673 OTUs. Further paired comparisons of 33 tumor tissues and normal adjacent tissues in terms of diversity metrics revealed significantly increased species richness (not adjusted two-tailed paired t test P value = 0.005) and related characteristics, Fisher’s diversity and the Chao1 index (P = 0.005), in normal samples. No significant differences (P > 0.05) in other indices of diversity were found between normal and tumor tissues. Further analysis revealed significant pairwise correlations (R = 0.73, P < 1E−5) of the species richness between tumors and paired normal tissues (Supplementary Fig. S2). A significant decrease in the number of species within the tumor was observed only in patients with a rich microbiome; no association was found when the tissue had low richness (Fig. 1A, Supplementary Fig. S3). There was also a trend (Fisher test P = 0.09) for females rather than males to have richer microbiomes in both tumor and normal tissues.Figure 1Associated changes in the richness and taxonomic structure of bacterial communities in tumor and matched normal tissues. (A) A significant difference in the number of identified spp. (richness) in tumors vs normal tissue was mostly observed in patients with a more diverse microbiome in normal tissue. The samples were sorted by richness and then divided into 3 equal richness groups (see dashed lines): low (L), medium (M), and high (H) number of species in the normal samples. A paired two-sided t test was used to evaluate whether tumor tissue was significantly different from normal tissue in each group. The p values produced from the test are shown in the figure. The sex of the patient is labeled as a bar on the Y axis with a rose (female) or blue (male) color. (B) Permutational multivariate analysis of variance of paired tumor-normal tissue samples in terms of bacterial species abundances revealed a significant association between the bacterial community structure and microbiome richness. Centroids of samples for each richness group are labeled H (High richness), M (Medium), and L (Low). (C) Permutational multivariate analysis of variance of paired tumor-normal samples from patients in the high-richness group revealed significant differences in the community structure between tumor and normal tissues. Centroids for tumor (red triangles) and normal (black circles) samples are labeled T and N, respectively. (D,E) Relative abundance of the most abundant phyla in normal and tumor tissues sorted by species richness in normal tissue. The phyla that are significantly more abundant in normal or tumor tissues are labeled with blue and red arrows, respectively. An increased abundance of Proteobacteriota and a decreased abundance of Actinobacteriota and Fusobacteriota were observed in tumors with low species richness. (F) Significant difference in the relative abundances of Bacteroidota, Fusobacteriota and Actinobacteriota between tumor and normal tissue from different richness groups. (G) Volcano plot of differentially abundant species identified by MaAsLin in normal and tumor tissues.The taxonomic structure of the bacterial community in normal and tumor samples also depended on species richness (Fig. 1B,C respectively). By applying permutational multivariate analysis of variance to all 33 paired tumor-normal samples, we observed significant differences in taxonomic composition when the paired samples were divided into 3 equal groups according to normal tissue richness: high (H), medium (M), and low (L) (Fig. 1D). The community structure of the H group differed from that of the M group and, especially, the L group at both the OTU level (Fig. 1D) and the Phylum level (Supplementary Fig. S4B). No difference was found between tumor and normal tissues when all the samples were analyzed together (Supplementary Fig. S4A). However, when examining each group separately, a significant difference in taxonomic structure between tumor and normal tissue was observed in the H-rich group of samples (Fig. 1E and Supplementary Fig. S4C) but not in the M or L group (Supplementary Fig. S4D). Due to inconsistent changes between the groups, more significant differences in the relative abundances of Phyla were observed if samples in the L and H richness groups were considered separately (Fig. 1F, Supplementary Figs. S5, S6). Normal and tumor tissues with low versus high species richness showed a significant increase in the abundance of Proteobacteriota (P = 0.001 and P = 0.0007) and a significant decrease in Fusobacteriota (P = 0.0003 and P = 0.002) (Fig. 1F). Tumors versus normals in H group (with high number of OTUs) showed increase of oral Bacteroidota and decrease of oral Actinobacteriota, and tumors versus normals in L group (low number of OTUs) showed decreased abundance of Bacteroidota (Fig. 1F) and increased abundance of Proteobacteriota (Supplementary Fig. S6).Associations between taxonomic structure and richness were also confirmed at the order level (Supplementary Fig. S7) and at the OTU level (Supplementary Fig. S8, Fig. 1G, Supplementary Fig. S9). Differentially abundant OTUs between normal and tumor tissues were found by comparison of 33 paired ACC tumor-normal tissues using MaAsLin2 (Fig. 1G) and LeFSe (Supplementary Figs. S8, S9, S10). Four putative species, including the typical inhabitants of the human oral cavity Veillonella parvula 45, Catonella morbi 46 and Alloscardovia omnicolens47, were abundant in normal tissue (Fig. 1G) and, specifically, in samples with a total number of species above the median value (Supplementary Fig. S9). The opposite relationship with species richness was found for 2 putative species differentially abundant in tumors (Supplementary Fig. S10): the gut bacterium Megamonas rupellensis48 and Bradyrhizobium centrosematis, which are known contaminants of purified and municipal water systems49. Both species were abundant in tumor tissue with low species richness. Moreover, Moraxella osloensis, which is implicated in ocular membrane inflammation50, was more abundant in tumor tissue but was not associated with species richness.Association of bacterial composition and abundances with overall patient survivalTo identify clinical factors and specific intratumoral species that might be associated with survival outcome, we analyzed tumor tissue from 50 patients (Supplementary Fig. S1). Analysis of clinical factors using the Cox proportional hazards model identified only one confounder, the ACC molecular subtype (ACC-I versus ACC-II), which was significantly associated with survival after multiple testing correction (log-rank Padj = 0.003) (Supplementary Table S2). Among the diversity characteristics, only the diversity inverse Simpson, diversity coverage, and Gini inequality indices showed statistical significance (log-rank Padj < 0.05), although they were not independent predictors of survival probability from the ACC molecular subtype (Supplementary Table S2).Among the topmost common and abundant species analyzed by the Cox proportional hazards model, only 1 OTU was significantly (P value = 8E−05) negatively associated with survival after adjustment of the log-rank p value for multiple comparisons (Fig. 2A). The OTU was classified taxonomically with 100% identity as a well-known commensal resident of the human gut51 Bacteroides thetaiotaomicron (B. theta). Importantly, the abundance of B. theta predicted survival probability independently of the ACC subtype, and the combination of B. theta abundance and ACC subtype increased the significance of the model from P = 2E−4 to P = 3E−5 according to the likelihood ratio test (Supplementary Table S3). Only 3 OTUs (Fig. 2A, Supplementary Table S2), taxonomically classified as typical inhabitants of the oral cavity (Granulicatella adiacens, Rothia mucilaginosa, and an unclassified species of Leptotrichia), had a significant positive association with survival (likelihood ratio test P, not adjusted; 0.03, 0.05, and 0.05, respectively). The abundances of the OTUs were significantly positively correlated with each other across the samples (rho ranged from 0.36 to 0.70), suggesting that the oral species can be a part of the same community. Each of the species had a significant negative correlation with B.theta abundance (rho ranged from − 0.27 to − 0.49), revealing their potential negative impact on the outgrowth of B. theta.Figure 2Grouping of ACC intratumoral microbiomes according to bacterial community type and abundance. (A) Bacterial species significantly associated with overall survival. (B) Grouping and ordering of ACCs according to the abundances of species associated with survival. The studied ACCs were subdivided into 2 groups: B. theta-positive and -negative. The former group was ordered by the sum of 3 OS-positive species abundances, from max to min, and the latter by B. theta abundance, from min to max. The B. theta-positive group of ACC tumors was dominated by the molecular subtype ACC-II, while the ACC-I subtype was significantly more common in B. theta-negative ACCs. (C) Association of B. theta− and B.theta + tumors with the age of the ACC patient. (D) Association of B. theta− and B.theta + tumors with ACC patient survival. (E) Supervised hierarchical clustering of the 54 most common spp. identified in ACCs reveals oral- and gut-type bacterial communities associated with B. theta-negative and -positive tumor groups, respectively. The oral-type bacterial community is enriched with common oral species, while the gut-type community is enriched with species that are common in the human gut. (F,G) Grouping of ACC intratumoral microbiomes according to oral and gut community abundances and associations of the groups according to richness; OH: high abundance of oral species, OL: low abundance of oral species, GH: high abundance of gut species, GL: low abundance of gut species. The OH and GL microbiomes had significantly greater richness than did the OL and GL microbiomes. (H) Permutational multivariate analysis of variance confirmed significant differences in the microbial composition of the proposed intratumoral microbiome groups (OH, OL, GH, and GL). (I) Significant association of the bacterial community taxonomic structure at the phylum level with the proposed grouping of the intratumoral microbiome in the ACC. Firmicutes, Actinobacteriota, and Fusobacteriota were more abundant in the B. theta-negative group, while the phyla Proteobacteria, Patescibacteria, and Cyanobacteria were more abundant in B. theta + tumors.To explore the associations of the species with other topmost common colonizers of the intratumoral microbiome, we grouped tumors and species in terms of B. theta and oral species abundances (Fig. 2B,C). We divided the tumors into B.theta + (positive) and B. theta− (negative) groups and then ordered the B. theta + tumors by B. theta abundance from the maximum to 0 (from right to left); and B. theta-negative samples were ordered by the sum of the 3 oral species abundances from min to max (Fig. 2B). Patients with B. theta− tumors were younger (Fig. 2C) and had better outcomes (Fig. 2D), and their tumors were more likely to be of the ACC-II myoepithelial molecular subtype (Fisher’s test P = 0.01). No associations were found between the B. theta− or B. theta + phenotype and primary tumor site (Pearson’s chi-square test P values of 0.23 and 0.37, respectively).Supervised hierarchical clustering was used to group the topmost common species according to similar abundance profiles across the ordered tumors (Fig. 2E). Similar profiles of a set of species across different conditions would suggest their existence as a community. The clustering revealed 2 major types of bacterial communities, referred to as oral type (O-type) and gut type (G-type). The O-type communities were dominated by species that are typical for the oral cavity, including Granulicatella adiacens, Rothia mucilaginosa, and Leptotrichia unclassified. The G-type communities were dominated by known common gut bacteria, including B. theta and Akkermansia muciniphila, which are well-known mucus layer degraders in the gut; putative species of Blautia and Lactobacillus gasserai; and other bacteria that are not typical for humans as hosts.Because the abundance of the O-type community varied in B. theta− tumors and decreased in tumors with low richness (Fig. 2F), we further subdivided the tumors into 2 subgroups (Fig. 2G): the Oral High (OH) group, which had a high abundance of oral species and high species richness in the intratumoral microbiome; and the Oral Low (OL) group, which had a low abundance of oral species and low species richness. In B. theta + tumors, the abundance of the G-type community also varied, but the abundance of gut species increased when the species richness was low (Fig. 2F). Therefore, the tumors were also subdivided into 2 subgroups: the gut high (GH) subgroup, which had a high abundance of gut species and low species richness in their intratumoral microbiomes, and the gut low (GL) subgroup, which had a low abundance of gut species and high species richness. The 4 identified subtypes of ACC intratumoral microbiomes, B. theta− OH and OL and B. theta + GH and GL, exhibited significant associations (nonparametric MANOVA P = 1e−04 by and PERMANOVA P = 0.001) with the taxonomic composition of the bacterial community (Fig. 2H,I). Species of Firmicutes, Actinobacteriota, and Fusobacteriota were more abundant in the B. theta-negative group (t test P = 0.06, P = 0.03, P = 0.10, respectively), while species of Proteobacteriota were significantly enriched in B. theta + tumors (t test P = 0.02), especially in samples characterized by low species richness (GH group), compared with the OH group (P = 0.002). Interestingly, Patescibacteria and Cyanobacteria, which are not common inhabitants of the oral cavity, had increased abundances in B. theta + tumors (t test P = 0.005 and P = 0.06, respectively).A significant positive correlation of species richness between paired tumor tissue and normal tissue in the present study (Supplementary Fig. S3) suggested that the dominance of certain oral or gut-associated species within the tumor might have occurred because bacteria from normal tissue colonized the developing tumor of the patient. Therefore, we explored the correlations of species abundances between tumor and normal tissues (Fig. 2E) in 33 paired samples to validate this hypothesis. Heatmaps of tumor and normal tissue paired with the same order of samples and OTUs as in Fig. 2E are presented in Fig. 3A,B, respectively. Correlation analysis for each species across 33 samples is provided in Fig. 3C. The results support the notion that tumors colonize oral species from normal tissue because their abundance in tumors significantly positively correlates with their abundance in normal tissue (mean R = 0.65, P < 0.02). The abundances of species composing the gut-type community in tumors were not correlated with those in normal tissues (mean R = 0.02; P > 0.1), suggesting de novo colonization of the tumor by non-oral, mainly gut-associated species.Figure 3Distinct characteristics of the gut- and oral-type communities in ACC tumors. (A,B) Heatmaps of the most common species in 33 pairs of ACC (A) and adjacent normal (B) samples. The species and samples in both heatmaps are shown in the same order as in Fig. 2C. Species of the gut-type community are abundant in ACC tumors but are weakly abundant or absent in their adjacent normal counterparts. (C) Pairwise Person correlation analysis showing a significant association between tumor and adjacent normal species in the oral-type community and no association for species in the gut-type community. (D,E) KEGG pathway and biological process enrichment in ACC cells significantly associated with the abundance of B. theta. (F) KEGG pathways significantly associated with the abundance of G. adiacens in ACC tumors.Biological processes correlated with the abundances of B. theta and G. adiacens in ACC tumorsTo identify biological processes correlated with the abundances of B. theta and G. adiacens in ACC tumors, we interrogated our previous RNA sequencing4 data for 44 ACC tumors via GSEA. The analysis (Fig. 3D) revealed only one KEGG pathway, biosynthesis of glycosphingolipids, that was significantly positively correlated with the B. theta abundance in ACC tumors (FDR = 0.03). Glycosphingolipids (GSLs) are known components of the cell membrane. They are composed of ceramide backbones and glycans that remain outside of the membrane and participate in cell adhesion, cell‒cell interactions, and migration52,53,54,55. The increased production of GSLs in tumors with abundant B. theta suggested that the availability of glycans for nutritional benefit might underlie the outgrowth of the bacterium. We also highlighted the activities of related processes, including assembly of the axoneme, cilium movement, and cilium-dependent cell motility (Fig. 3E)56, among the top biological processes positively correlated with B. theta. In contrast, activities associated with the normal structure and organization of the salivary gland, such as gland morphogenesis and basement membrane organization, were negatively correlated with B. theta abundance (Fig. 3E).Several molecular pathways known to be significantly more active in the ACC-II molecular subtype4, including extracellular matrix receptor interactions, focal adhesion, Hedgehog, and TGF-beta signaling pathways, were positively correlated with G. adiacens abundance (Fig. 3F). Likewise, processes associated with cell proliferation, a known signature of the ACC-I molecular subtype, were negatively correlated with G. adiacens abundance.Comparison of intratumoral microbiomes in ACC patients with oral and head and neck cancersThe bacterial communities identified in ACC may resemble those found in oral squamous cell carcinoma (OSCC) and head and neck squamous cell carcinoma (HNSC), as all three cancers originate in or near the oral cavity. To identify the bacterial genera that cancers may share with ACC, we examined previously published datasets of the intratumoral microbiomes of these cancers. We selected the most common genera from each of the three OTU tables and clustered each table using unsupervised hierarchical clustering with the same parameters (Fig. 4A–C). Clustering identified 2 communities of bacterial Genera (Community 1 and Community 2) associated with each cancer. Many Genera were shared between ACC and OSCC (36%) and between ACC and HNSC (34%) and clustered in Community 1, which is also referred to as shared. Importantly, 6 Genera of the healthy oral microbiome, Neisseria, Leptotrichia, Actinomyces, Streptococcus, Rothia, and Veillonella, were shared among all three cancers and clustered close to each other in shared Community 1. Like that in ACC, the bacterial community in OSCC and HNSC tumors was enriched and characterized by a more diverse intratumoral microbiome (Fig. 4D,F, respectively), a less aggressive phenotype (Fig. 4B,C, respectively), and better overall patient survival (Fig. 4E,G, respectively). Thus, the shared oral genera may be considered biomarkers of healthier tumor microbiomes in all three cancers. In contrast, Community 2 was composed of distinct bacteria from each of the three cancers and was associated with a more aggressive phenotype and worse patient survival.Figure 4Comparative analysis of the intratumoral microbiomes of ACC patients with those of OSCC and HNSC patients. (A–C) Unsupervised hierarchical clustering of tumors in the ACC (A), OSCC (B), and HNSC (C) cohorts identified 2 bacterial communities (rows) in each cohort: Community 1 populated by shared oral genera (names are colored in red) and Community 2 populated by genera that are mainly specific for each cohort (names are colored in black). The names of the genera identified in all three cohorts are given in red bold. The columns of each heatmap denote tumors that are clustered together into major tumor groups (top of each heatmap) with similar abundance profiles of the intratumoral communities. The bars at the bottom of each heatmap show characteristics associated with the tumor groups. In the ACC (A), clustering revealed 2 tumor groups that recapitulated the supervised clustering (Fig. 2B,E) based on intratumoral B. theta abundance. Group 1 ACC tumors had many healthy oral genera but less abundant Bacteroides (t-test P = 0.0006) and other gut-associated Genera. The aggressive ACC I molecular subtype was also less common in patients in the other groups (Fisher’s test P = 0.14). In the OSCC cohort (B), Group 1 tumors had less abundant Fusobacterium (t test P = 0.0001), a significantly reduced frequency of TP53 mutations (Fisher’s test P = 0.037) and disease recurrence (Fisher’s test P = 0.044). Bacteroides and Fusobacterium are labeled by red rectangles in ACC and OSCC, respectively, as potential biomarkers of a more aggressive phenotype. In HNSC (C), clustering revealed 3 tumor groups (Group 1, Group 2, and Group 3). Group 1 HNSC tumors were less likely than Group 3 HNSC tumors to be at an advanced stage (Fisher’s test P = 0.01) or to spread into the lymph nodes (Fisher’s test P = 0.05). No difference was found between Group 1 and Group 2 tumors. (D,E) Community 1 in OSCC patients was associated with a more diverse microbiome (D) and better overall patient survival (E). (F,G) Community 1 in HNSC was also associated with a more diverse microbiome (F) and better overall patient survival (G).
