Clinical characteristics of the cohorts
Fresh 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 (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 study [4]. 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%).
The taxonomic structure of the bacterial community in ACC tumors and adjacent normal tissue is associated with bacterial richness
16S 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 (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, 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.
The taxonomic structure of the bacterial community in tumor and normal samples also depended on species richness. 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. 1B). 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. 1B) and the Phylum level (Fig. S4B). No difference was found between tumor and normal tissues when all the samples were analyzed together (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. 1C and Fig. S4C) but not in the M or L group (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, Fig. 1G, Fig. S5, and Fig. 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, Fig. 1F) 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 (Fig. 1G), and tumors versus normals in L group (low number of OTUs) showed decreased abundance of Bacteroidota (Fig. 1G) and increased abundance of Proteobacteriota (Fig. S6).
Associations between taxonomic structure and richness were also confirmed at the order level (Fig. S7) and at the OTU level (Fig. S8, Fig. 1J, Fig. S9). Differentially abundant OTUs between normal and tumor tissues were found by comparison of 33 paired ACC tumor-normal tissues using MaAsLin2 (Fig. 1J) and LeFSe (Fig. S8, Fig. S9, Fig. S10). Four putative species, including the typical inhabitants of the human oral cavity Veillonella parvula [45], Catonella morbi [46] and Alloscardovia omnicolens [47], were abundant in normal tissue (Fig. 1J) and, specifically, in samples with a total number of species above the median value (Figure S9). The opposite relationship with species richness was found for 2 putative species differentially abundant in tumors (Fig S10): the gut bacterium Megamonas rupellensis [48] and Bradyrhizobium centrosematis, which are known contaminants of purified and municipal water systems [49]. Both species were abundant in tumor tissue with low species richness. Moreover, Moraxella osloensis, which is implicated in ocular membrane inflammation [50], was more abundant in tumor tissue but was not associated with species richness.
Association of bacterial composition and abundances with overall patient survival
To identify clinical factors and specific intratumoral species that might be associated with survival outcome, we analyzed tumor tissue from 50 patients (Figure 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) (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 (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 gut [51] 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 (Table S3). Only 3 OTUs (Fig. 2A, 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.
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 and 2C). 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, Fig. 2I). 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 (Figure 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 and Fig. 3B, 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.
Biological processes correlated with the abundances of B.theta and G. adiacens in ACC tumors
To identify biological processes correlated with the abundances of B. theta and G. adiacens in ACC tumors, we interrogated our previous RNA sequencing [4] 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 migration. [52–55]. The increased production of GSLs in tumors with abundant B.theta suggested that the availability of glycans for nutritional benefit may 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 subtype [4], 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 cancers
The 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, Fig. 4B, Fig. 4C). 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 (Figs. 4D and 4F, respectively), a less aggressive phenotype (Figs. 4B and 4C, respectively), and better overall patient survival (Figs. 4E and 4G, 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.