Abstract
Esophageal adenocarcinoma (EAC) is rising in incidence, and established risk factors do not explain this trend. Esophageal microbiome alterations have been associated with Barrett's esophagus (BE) and dysplasia and EAC. The oral microbiome is tightly linked to the esophageal microbiome; this study aimed to identify salivary microbiome-related factors associated with BE, dysplasia, and EAC.
Clinical data and oral health history were collected from patients with and without BE. The salivary microbiome was characterized, assessing differential relative abundance of taxa by 16S rRNA gene sequencing and associations between microbiome composition and clinical features. Microbiome metabolic modeling was used to predict metabolite production.
A total of 244 patients (125 non-BE and 119 BE) were analyzed. Patients with high-grade dysplasia (HGD)/EAC had a significantly higher prevalence of tooth loss (P = 0.001). There were significant shifts with increased dysbiosis associated with HGD/EAC, independent of tooth loss, with the largest shifts within the genus Streptococcus. Modeling predicted significant shifts in the microbiome metabolic capacities, including increases in L-lactic acid and decreases in butyric acid and L-tryptophan production in HGD/EAC.
Marked dysbiosis in the salivary microbiome is associated with HGD and EAC, with notable increases within the genus Streptococcus and accompanying changes in predicted metabolite production. Further work is warranted to identify the biological significance of these alterations and to validate metabolic shifts.
There is an association between oral dysbiosis and HGD/EAC. Further work is needed to establish the diagnostic, predictive, and causal potential of this relationship.
Introduction
Esophageal adenocarcinoma (EAC) has seen a dramatic rise in incidence in the past several decades, is often diagnosed at advanced stages, and is associated with poor survival (1, 2). The factors that drive EAC remain incompletely understood. Barrett's esophagus (BE) is the precursor lesion to EAC, but the overwhelming majority of BE patients do not progress to EAC. Established EAC risk factors, including gastroesophageal reflux disease (GERD) and obesity, do not fully explain the rise in its incidence (3, 4).
Increasing evidence suggests that the microbiome plays an important role in modifying the risk of a variety of epithelial cancers (5–8) as well as in modulating the response to treatment (9–11). Changes in the esophageal microbiome have been observed in EAC and with progression from BE to EAC (12, 13), raising the possibility that bacteria contribute to esophageal neoplasia. Reliable sampling of the esophageal microbiome, however, requires invasive procedures. A more accessible “window” to the esophageal ecosystem is the oral microbiome, which correlates with it (14). A small study of the tumor-associated microbiome in EAC found a high prevalence of domination by oral flora such as Streptococcus (12), pointing to a link between the oral microbiome and EAC. Oral microbiome alterations have been associated with future risk of EAC (15), and differences in the oral microbiome of BE patients were described previously in a small study of 49 patients (16). Alterations in the oral microbiome have also been associated with poor oral health (17), which was in itself associated with an increased risk of EAC in a recent analysis (18). It remains unclear how oral dysbiosis and poor oral health interact in their association with EAC. Finally, little is known with regard to oral microbiome alterations associated with dysplasia and EAC in BE patients. A clearer understanding of these oral microbiome changes could identify factors that may drive the progression of neoplasia, representing novel therapeutic targets.
Here, we profiled the salivary microbiome from 250 patients with various stages of BE and EAC who were undergoing upper endoscopy. We identify multiple characteristics of the oral microbiome associated with stages of dysplasia and early cancer in BE and show that they are independent of oral health. Using metabolic modeling, we predict metabolite profiles associated with alterations in BE, suggesting a mechanistic role for microbially produced metabolites. Finally, we show that the salivary microbiome offers a mild improvement in diagnostic accuracy compared with models based on clinical risk factors to identify patients with high-grade dysplasia (HGD) or EAC. Our results demonstrate the potential of studying the oral microbiome in the context of progression to EAC.
Materials and Methods
Study design
A total of 250 patients with and without BE undergoing upper endoscopy at Columbia University Irving Medical Center (New York, NY) were prospectively enrolled from February 2018 through February 2019. Patients were ≥18 years old and scheduled to undergo endoscopy for clinical indications. Patients were excluded if they had a concurrently scheduled colonoscopy, had a history of gastric or esophageal surgery, a history of esophageal squamous cell cancer, or use of antibiotics, steroids, or other immunosuppressants in the 3 months prior to the procedure. This study was approved by the Columbia University Institutional Review Board. All patients provided written informed consent. The study was carried out in accordance with the principles of the Declaration of Helsinki.
Data were collected on patient demographics and anthropometrics (to calculate body mass index, BMI) as well as clinical information including medical history, history of gastroesophageal reflux disease (GERD; defined as experiencing frequent heartburn or fluid regurgitation), medication use at time of enrollment (with specific notation of daily use of proton pump inhibitors (PPI), histamine-2 receptor antagonists, statins, and daily use of aspirin and nonsteroidal anti-inflammatory drugs), alcohol history, and smoking history (ever smoking defined as having smoked >100 lifetime cigarettes). Data were collected on self-reported oral health and hygiene. Tooth loss was assessed using categories adapted from Borningen and colleagues (17): all or most of natural adult teeth, partial plates or implants, full upper dentures or implants, full lower dentures or implants, full upper and lower dentures or implants. Data were also collected on tooth brushing and mouthwash use.
Patients did not eat or drink after midnight prior to the endoscopy and saliva collection; saliva was collected prior to the endoscopy. Patients were categorized as BE if they had a history of endoscopically suspected BE with intestinal metaplasia on esophageal biopsies. BE patients were further categorized based on the highest degree of neoplasia ever [no dysplasia (NDBE), indefinite for dysplasia (IND), low-grade dysplasia (LGD), high-grade dysplasia (HGD), and adenocarcinoma (EAC)].
Microbiome sequencing and analysis
DNA was extracted using Qiagen MagAttract PowerMicrobiome DNA/RNA. A positive control (ZymoBIOMICS Microbial Community Standard, D6300) and two negative controls (sterile water) were included in each extraction plate. The 16S rRNA V3-V4 region was amplified using Illumina adapter-ligated primers (19). The Illumina Nextera XT v2 index sets A–D were used to barcode sequencing libraries. Libraries were sequenced on an Illumina MiSeq using the v3 reagent kit (600 cycles) and a loading concentration of 12 pmol/L with 10% PhiX spike-in. Sequences were assigned to operational taxonomic units (OTU) using USEARCH (biorxiv 2016.09.09.074161) with ≥ 97% sequence homology. Taxonomic assignments for the OTUs were based on the Human Oral Microbiome Database (20). Any subsequently unassigned OTUs were assigned by referencing the Ribosomal Database Project (21). Samples were subsampled to 10,000 reads to compare across even sequencing depths while minimizing data loss. Analysis of negative controls did not reveal significant contamination. Five patients were excluded after sequencing because of relatively low sequencing depth with <10,000 total reads per sample. The median read count for the full cohort was >33,000. One patient was excluded because of a history of both EAC and esophageal squamous cell carcinoma.
Microbiome metabolic modeling of oral microbial communities
Microbiome metabolic modeling was performed using the Microbiome Modeling Toolbox (COBRA toolbox commit: 71c117305231f77a0292856e292b95ab32040711; refs. 22, 23) and the AGORA metabolic models (AGORA 1.02) (24). All computations were performed in MATLAB version 2019a (Mathworks, Inc.), using the IBM CPLEX (IBM, Inc.) solver. We first matched species detected by our microbial sequencing analysis with the ones present in AGORA (24). Because AGORA metabolic models are available at the strain level, we generated species-level models using the createPanModels.m function of the Microbiome Modeling Toolbox (MMT) (22) as previously described (25). To increase the number of species represented in our microbiome models, we chose genus-level representative models for abundant microbes present in the oral cavity with >5% relative abundance in more than 10 samples. There were six species without a corresponding metabolic model, and these were either grouped with similar species or excluded from the analyses (see Supplementary Table S1 for details).
We then used the mgPipe.m automated pipeline of the MMT to build and interrogate sample-specific microbiome metabolic models. Briefly, for each sample, personalized microbiome models are created by joining species-level metabolic models using the compartmentalization technique (26); a lumen compartment enabling microbial metabolic interactions is added, as well as additional input and output compartments, allowing microbiome intake and secretion of metabolites. Altogether, our microbiome models included 160 microbial species with an average of 50 species for each sample and a maximum of 69. As constraint-based metabolic modeling benefits from a specification of the metabolic environment such as media and carbon source availability (27), we applied a “western diet” (28) to each sample in the form of constraints on the metabolites uptake reactions (28). Finally, to obtain metabolic predictions, we used the Net Maximal Production Capabilities (NMPC) through the mgPipe pipeline (22) to provide predictions of the metabolite secretion profile of each sample. To detect significant changes in NMPC distributions between cases and controls a Mann–Whitney U test was performed for each retained NMPCs. Only NMPCs which were present in at least 10% of the cases and had at least a value of 0.01 were retained for the significance analysis. FDR correction using the Benjamini–Hochberg procedure was applied.
Statistical analysis
The primary groups of comparison were BE patients with HGD/EAC, nondysplastic BE (NDBE) and non-BE controls. Grouping high-grade dysplasia and intramucosal adenocarcinoma together as advanced neoplasia reflects common practice as well as clinical guidelines for treatment (29). There is extremely low interobserver agreement (even among expert gastrointestinal pathologists) for the diagnosis of LGD (30–32), as inflammation-induced cytologic atypia mimics the findings of LGD. As a result, although estimates of cancer risk for LGD are relatively low on average (31), these estimates vary widely, thus making interpretations of findings for this group challenging. Patients with LGD or indefinite for dysplasia were included in analyses assessing for alterations in the oral microbiome across the entire BE neoplastic spectrum. Although there was no a priori reason to suspect that endoscopic therapy would have altered the salivary microbiome, comparisons were made between those patients with LGD or worse who had (n = 78) and had not (n = 10) received prior endoscopic therapy. There were no differences in alpha diversity (P = 0.16), no evidence of clustering on beta diversity analyses (ANOSIM; P = 0.13), and no differentially abundant taxa. Thus, treated and untreated patients were grouped together for all analyses.
Categorical variables were compared across groups using Fisher exact tests. Continuous variables were analyzed using t tests or rank sum tests as appropriate, with ANOVA and Kruskal–Wallis tests for ≥3 groups. For purposes of analyses, tooth loss was dichotomized as having all or most of natural adult teeth (yes/no). As indicated above, all subjects who did not have all or most of natural adult teeth had dental hardware; thus, analyses could not distinguish between having all or most of adult natural teeth and the presence or absence of dental hardware. Multivariable logistic regression was performed to assess the association between tooth loss and advanced neoplasia, adjusted for known EAC risk factors (age, sex, GERD, BMI, and smoking).
Alpha diversity was evaluated using the Shannon diversity index and beta diversity using weighted UniFrac (33) distances. Groups were compared using both permutational multivariate analysis of variance (PERMANOVA) for predicted metabolite profiles and analysis of similarities (ANOSIM) for microbial compositions. To find differential abundances between study groups, the ALDEx2 (34) R package was used. For differential abundance analyses, only OTUs present in at least 5% of all samples were included to allow for more meaningful comparisons. ALDEx2 was used to compare worst histologic grades of BE as an ordinal variable in a generalized linear model and to assess correlation of BE-associated OTUs with neoplastic progression using aldex.corr to treat worst histologic grade as a continuous variable. ALDEx2 was also used to find significance for differentially abundant taxa in a multivariate model with both advanced neoplasia and tooth loss.
Generalized linear models were used to assess the differential relative abundance of bacterial taxa in advanced neoplasia, adjusted for tooth loss. Multivariable logistic regression was performed to detect associations between advanced neoplasia and microbiome composition (represented by its top five principal coordinates), adjusted for EAC risk factors (age, sex, race, BMI, smoking, GERD). Supervised machine learning was used to classify patients with advanced neoplasia using the LightGBM package (35). Three models were created: (i) EAC risk factors alone (age, sex, race, BMI, smoking, GERD); (ii) microbiome features alone; and (iii) EAC risk factors and microbiome features together. Model parameters were optimized per fold in 10-fold cross-validation, with strict train–test sterility. The output of the models was predicted probabilities of whether a patient has advanced neoplasia or no BE, with the goal of identifying the patients at the highest risk of mortality from EAC.
Post hoc power calculations demonstrated that, assuming alpha = 0.05, the study sample size (119 BE subjects, 125 controls) had 87% power to detect at 0.4 SD difference in Shannon index between BE and controls. Considering the 200 most abundant taxa, the sample size had 80% and 90% power to detect statistically significant (Bonferonni-adjusted P < 0.1) differences in relative abundance between the groups in 18 and 13 taxa, respectively. All statistical analyses were performed in Python or R. Statistical significance was defined as P < 0.05. Differential abundance analyses were corrected for multiple comparisons using the Benjamini–Hochberg procedure, and corrected statistical significance was defined as P < 0.1. 95% confidence intervals for areas under the receiver operating characteristic curve (AUROC) were calculated using the DeLong method using pROC (36).
Data availability
16S rRNA gene sequencing files were uploaded to NCBI Sequence Read Archive (PRJNA785879).
Results
Oral microbial composition from a large endoscopy cohort
We recruited 250 adult patients undergoing upper endoscopy and characterized their oral microbiome using 16S rRNA gene sequencing. Six subjects were excluded from analyses: one also had a history of esophageal squamous cell cancer in addition to EAC, and 5 had <10,000 reads across all OTUs. A total of 244 patients were included in the analyses: 125 controls without BE, and 119 BE patients [20 with nondysplastic BE, 11 indefinite for dysplasia, 10 LGD, 54 HGD, and 24 intramucosal (T1a) adenocarcinoma]. The indications for endoscopy in the non-BE controls are shown in Supplementary Table S2. Patients with BE were more likely to be older (t test P < 0.001), male (Fisher exact P < 0.001), white (P = 0.001), or ever-smokers (defined as ≥100 lifetime cigarettes smoked) (P = 0.003). They were also more likely to have GERD (P < 0.001), to be treated with PPI (P < 0.001), to take aspirin (P < 0.001), and to have a higher BMI (P < 0.001) (Table 1; Supplementary Table S3). There was no significant difference in the use of mouthwash between BE and non-BE patients (P = 0.61), but non-BE patients were more likely to brush their teeth at least daily (98% vs. 92%, P = 0.03). Compared with non-BE, a significantly higher proportion of patients with BE had tooth loss (63% vs. 42%, P = 0.001), largely due to an increase in tooth loss in patients with HGD/EAC (66%) compared with nondysplastic BE (40%) and non-BE (42%; Fisher exact P = 0.001) (Fig. 1). Older age (per year, adjusted OR 1.06; 95% CI, 1.04–1.08) and a history of smoking (adjusted OR 2.12; 95% CI, 1.07–4.19) were independently associated with tooth loss (Supplementary Table S4). In multivariable analyses adjusting for EAC risk factors (age, male sex, white race, and GERD), tooth loss was associated with a nonsignificant increased risk of HGD/EAC (vs. non-BE, adjusted OR 1.49; 95% CI, 0.96–2.47) (Supplementary Table S5). Established EAC risk factors were independently associated with HGD/EAC even after controlling for daily tooth brushing, use of mouthwash, and presence of tooth loss, whereas these measures of oral health and hygiene were not independently associated with HGD/EAC (P = 0.21, 0.57, and 0.34, respectively).
Patient characteristics.
. | All (n = 244) . | Non-BE (n = 125) . | BE (n = 119) . | P value . |
---|---|---|---|---|
Age, years; mean (SD) | 57.8 (18.7) | 50.9 (18.7) | 65.0 (15.8) | <0.001 |
Sex, male; N (%) | 140 (57%) | 45 (36%) | 95 (80%) | <0.001 |
Ever-smokers; N (%) | 103 (42%) | 41 (33%) | 62 (52%) | 0.003 |
BMI mean (SD) | 27.5 (6.7) | 27.0 (7.9) | 28.1 (5.1) | <0.001 |
Race, white; N (%) | 222 (90%) | 105 (84%) | 117 (98%) | 0.001 |
GERD; N (%) | 167 (68%) | 64 (51%) | 103 (87%) | <0.001 |
PPI use; N (%) | 148 (60%) | 41 (33%) | 107 (90%) | <0.001 |
Aspirin use; N (%) | 77 (31%) | 25 (20%) | 52 (44%) | <0.001 |
Oral health and hygiene | ||||
Tooth loss; N (%) | 127 (52%) | 52 (42%) | 75 (63%) | <0.001 |
Tooth brushing frequency, ≥daily; N (%) | 233 (95%) | 123 (98%) | 110 (92%) | 0.03 |
Mouthwash frequency, ≥daily; N (%) | 139 (56%) | 69 (55%) | 70 (59%) | 0.61 |
. | All (n = 244) . | Non-BE (n = 125) . | BE (n = 119) . | P value . |
---|---|---|---|---|
Age, years; mean (SD) | 57.8 (18.7) | 50.9 (18.7) | 65.0 (15.8) | <0.001 |
Sex, male; N (%) | 140 (57%) | 45 (36%) | 95 (80%) | <0.001 |
Ever-smokers; N (%) | 103 (42%) | 41 (33%) | 62 (52%) | 0.003 |
BMI mean (SD) | 27.5 (6.7) | 27.0 (7.9) | 28.1 (5.1) | <0.001 |
Race, white; N (%) | 222 (90%) | 105 (84%) | 117 (98%) | 0.001 |
GERD; N (%) | 167 (68%) | 64 (51%) | 103 (87%) | <0.001 |
PPI use; N (%) | 148 (60%) | 41 (33%) | 107 (90%) | <0.001 |
Aspirin use; N (%) | 77 (31%) | 25 (20%) | 52 (44%) | <0.001 |
Oral health and hygiene | ||||
Tooth loss; N (%) | 127 (52%) | 52 (42%) | 75 (63%) | <0.001 |
Tooth brushing frequency, ≥daily; N (%) | 233 (95%) | 123 (98%) | 110 (92%) | 0.03 |
Mouthwash frequency, ≥daily; N (%) | 139 (56%) | 69 (55%) | 70 (59%) | 0.61 |
P = Wilcoxon rank sum or Fisher exact P for the difference between BE and non-BE.
Tooth loss is significantly more common in advanced neoplasia. A significantly higher proportion of patients with advanced neoplasia (high-grade dysplasia or esophageal adenocarcinoma; n = 78) had tooth loss as compared with non-BE (n = 125) and nondysplastic BE (n = 20) patients combined (Fisher exact P = 0.001). P = Fisher exact.
Tooth loss is significantly more common in advanced neoplasia. A significantly higher proportion of patients with advanced neoplasia (high-grade dysplasia or esophageal adenocarcinoma; n = 78) had tooth loss as compared with non-BE (n = 125) and nondysplastic BE (n = 20) patients combined (Fisher exact P = 0.001). P = Fisher exact.
The oral microbiome of BE patients is progressively altered with dysplastic changes
To assess whether the salivary microbiome is associated with stages of neoplastic progression, we focused our analyses on comparisons between three groups: non-BE (n = 125), nondysplastic BE (n = 20), and HGD/EAC (n = 78). Across these stages, we found significantly different alpha diversity (Shannon: Kruskal–Wallis P = 0.005; Fig. 2A; Simpson: P = 0.0029; Supplementary Fig. S1). Compared with patients without BE, the alterations in alpha diversity were more pronounced in patients with HGD/EAC than in those with nondysplastic BE (non-BE vs. nondysplastic BE, Mann–Whitney P = 0.11; non-BE vs. advanced neoplasia, P = 0.0006). There was no significant difference in alpha diversity comparing nondysplastic BE with HGD/EAC (P = 0.23). We further found that the oral microbiome from patients with HGD/EAC tended to cluster separately than the rest of the cohort (weighted UniFrac, ANOSIM P < 0.001; Fig. 2B). Similar results were found when including all the subjects in their individual groups (non-BE, nondysplastic BE, IND, LGD, HGD, and EAC; Supplementary Fig. S2). Our results indicate that the salivary microbiome alterations observed in HGD/EAC are reflected both in the diversity of each individual's microbiome (alpha diversity) and in compositional differences between individuals (beta diversity).
A microbial signature of BE. A, Patients with advanced neoplasia have significantly reduced alpha diversity compared with non-BE patients. Kruskal–Wallis overall P = 0.005. P – Mann–Whitney U test. B, Weighted UniFrac PCoA demonstrated significant clustering of patients with advanced neoplasia (ANOSIM P = 0.001). NDBE, nondysplastic BE; advanced neoplasia, high-grade dysplasia or esophageal adenocarcinoma; ellipse, 2 standard deviation sigma ellipse.
A microbial signature of BE. A, Patients with advanced neoplasia have significantly reduced alpha diversity compared with non-BE patients. Kruskal–Wallis overall P = 0.005. P – Mann–Whitney U test. B, Weighted UniFrac PCoA demonstrated significant clustering of patients with advanced neoplasia (ANOSIM P = 0.001). NDBE, nondysplastic BE; advanced neoplasia, high-grade dysplasia or esophageal adenocarcinoma; ellipse, 2 standard deviation sigma ellipse.
We next checked whether specific microbes were associated with BE stages of dysplasia and EAC. To do this, we first compared all BE and non-BE to establish the key BE-associated OTUs using ALDEx2 (34). A total of 26 OTUs were identified as differentially abundant (P < 0.05, FDR corrected at 0.1; Fig. 3). To assess whether the dysbiotic signature associated with BE is more pronounced with dysplastic changes, we then analyzed these 26 OTUs in different stages of BE, from no dysplasia to EAC. There was a significant association (P < 0.05, FDR corrected at 0.1) for 23 of the 26 taxa, with a notable shift in composition in the transition from LGD to HGD (Fig. 3). This transition in composition from LGD to HGD is consistent with our prior observations of esophageal microbiome alterations in BE (13). The taxonomic alterations were notable for the increased relative abundance of several Streptococcus species. Streptococci form biofilms in the oral cavity (37–39).
Increased oral dysbiosis with progressive dysplasia. Shifts in the oral microbiome compared with non-BE patients (shown as ALDEx2 effect sizes) were more pronounced with progression from no dysplasia to EAC, particularly notable in patients with high-grade dysplasia and EAC. Bolded OTUs were significantly associated with neoplastic progression (P < 0.05, FDR corrected at 0.1) NDBE, nondysplastic Barrett's esophagus (n = 20); Indef, indefinite for dysplasia (n = 11); LGD, low-grade dysplasia (n = 10); HGD, high-grade dysplasia (n = 54); EAC, esophageal adenocarcinoma (n = 24).
Increased oral dysbiosis with progressive dysplasia. Shifts in the oral microbiome compared with non-BE patients (shown as ALDEx2 effect sizes) were more pronounced with progression from no dysplasia to EAC, particularly notable in patients with high-grade dysplasia and EAC. Bolded OTUs were significantly associated with neoplastic progression (P < 0.05, FDR corrected at 0.1) NDBE, nondysplastic Barrett's esophagus (n = 20); Indef, indefinite for dysplasia (n = 11); LGD, low-grade dysplasia (n = 10); HGD, high-grade dysplasia (n = 54); EAC, esophageal adenocarcinoma (n = 24).
We assessed whether smoking history was associated with salivary microbiome composition. Only a small number of subjects in the study were current smokers (7/244, 2.9%), precluding meaningful analyses of active smoking and the salivary microbiome. Comparing ever- versus never-smokers, there was no significant difference in alpha diversity (Shannon index, P = 0.08), evidence of clustering on beta diversity analyses (ANOSIM P = 0.08), and no differentially abundant taxa.
We also assessed whether individuals in the control group who had celiac disease may have been driving some of the observed differences, as some prior studies have reported oral microbiome alterations in celiac disease patients (40). In comparisons of celiac disease versus nonceliac controls, we found no differences in composition based on alpha diversity (Shannon index, P = 0.12) and beta diversity analyses (ANOSIM P = 0.87), and there were no significantly differentially abundant taxa. We also repeated the analyses comparing HGD/EAC to controls, excluding those with celiac disease, and there were no substantive changes in the results. HGD/EAC patients again had significantly reduced alpha diversity (Shannon index, P = 0.009), and there was significant clustering on beta diversity analyses (ANOSIM P < 0.001). There were 8 significantly differentially abundant taxa, all of which were also identified as differentially abundant in the original comparisons with the full control group.
The salivary microbiome is associated with HGD and EAC even when controlling for tooth loss
Tooth loss is known to be strongly associated with oral microbiome composition (17), and there was an increased proportion of HGD/EAC patients with tooth loss (Fig. 1). Comparing patients who did and did not have all or most of their natural adult teeth, we found that those with tooth loss had lower alpha diversity (Shannon, Mann–Whitney P = 0.001; Supplementary Fig. S3A), and that oral microbiomes from both groups clustered separately (weighted UniFrac, ANOSIM P = 0.003; Supplementary Fig. S3B). We also identified 29 OTUs that had significantly different abundance between patients with and without tooth loss (ALDEx2 P < 0.05, FDR corrected at 0.1; Fig. 4). As poor oral health is associated with oral dysbiosis, this raised the question of whether the oral microbiome is associated with HGD/EAC independent of tooth loss.
Oral microbes are independently associated with advanced neoplasia and tooth loss. Differentially abundant taxa in patients with advanced neoplasia (high-grade dysplasia or esophageal adenocarcinoma) compared with non-BE (left) and in patients with and without tooth loss (right). Many taxa were associated with both neoplasia and tooth loss, yet most of these remained significantly differentially abundant after adjusting for tooth loss and advanced neoplasia, respectively (ALDEx2 P < 0.05, FDR corrected at 0.1; denoted by asterisks).
Oral microbes are independently associated with advanced neoplasia and tooth loss. Differentially abundant taxa in patients with advanced neoplasia (high-grade dysplasia or esophageal adenocarcinoma) compared with non-BE (left) and in patients with and without tooth loss (right). Many taxa were associated with both neoplasia and tooth loss, yet most of these remained significantly differentially abundant after adjusting for tooth loss and advanced neoplasia, respectively (ALDEx2 P < 0.05, FDR corrected at 0.1; denoted by asterisks).
We first examined whether salivary microbiome composition as a whole is associated with HGD/EAC independently of tooth loss. We therefore calculated microbiome principal coordinates (PCo) using weighted UniFrac distances and used the top five PCos, which represented two thirds of the variance in microbiome composition. We then used multivariable logistic regression and found that PCo2 (explaining 22% of microbiome variation) and PCo4 (4.6%) were independently associated with HGD/EAC (P < 0.001 and P = 0.004, respectively). We then added the major EAC risk factors (age, sex, race, BMI, GERD history, and smoking history) to the model, and found that PCo2 remained independently associated with advanced neoplasia (P = 0.004), suggesting that salivary microbiome composition represents a potential novel independent risk factor for EAC. Adding tooth loss to the model did not alter the association between PCo2 and HGD/EAC (P = 0.004), and in this model, tooth loss was not independently associated with HGD/EAC (P = 0.12). To address the possibility that gastric acid–suppressing medications may act as a confounder, we added PPI use to the model, and PCo2 remained significantly associated with HGD/EAC (P = 0.03). Our results suggest that the association of tooth loss with advanced neoplasia is mediated through the oral microbiome.
We next assessed whether associations between specific oral taxa and HGD/EAC are independent of tooth loss. Of the 33 taxa associated with HGD/EAC (N = 78) versus non-BE (N = 125; ALDEx2 P < 0.05; FDR corrected at 0.1), 18 were also associated with tooth loss (Fig. 4). After adjusting for tooth loss in a generalized linear model, 20 of these taxa remained significantly associated with advanced neoplasia (ALDEx2 P < 0.05, FDR corrected at 0.1). Notably, the four OTUs with the greatest increase in relative abundance in HGD/EAC were all assigned to the genus Streptococcus, and the increased abundance of these Streptococcus OTUs in HGD/EAC was independent of tooth loss. This corresponds with previous studies that found that the tumor-associated microbiome in EAC is often dominated by Streptococcus species (12).
Metabolic modeling predicts distinct metabolic secretion capabilities in HGD and EAC
Metabolite production by microbial communities is an important modality by which the microbiome affects the host. In order to assess if microbially produced metabolites might play a role as a driver or biomarker of neoplasia, we used microbiome community-scale metabolic models to predict metabolite secretion by the microbiome for every sample. We found significant clustering of predicted metabolite profiles comparing HGD/EAC cases with non-BE subjects (PERMANOVA P = 0.001). Using principal components analysis, we found a notable shift in the second component (15% explained variance; Mann–Whitney P = 0.0003; Fig. 5A) Forty-four predicted metabolites had significantly altered abundance (P < 0.05, FDR corrected at 0.1) in HGD/EAC (Fig. 5B). Notable alterations included increased predicted levels of L-lactic acid (P = 0.023), a by-product of aerobic glycolysis, a hallmark of cancer that can contribute to neoplasticity (41); and 2-ketobutyric acid (P = 0.033), previously reported to support mitochondrial respiration and cell proliferation (42). We also predicted that HGD/EAC features a decrease in butyric acid (P = 0.0089), a key promoter of gut homeostasis that was previously shown to be depleted in colon cancer and inflammatory bowel disease (43–45); and a decrease in L-tryptophan (P = 0.0017) (Fig. 5C). Circulating levels of tryptophan have been inversely associated with colon cancer risk (46), and melatonin, a by-product of L-tryptophan metabolism, is under investigation for EAC prevention (47). The pattern of the shifts we predict is therefore consistent with the potential promotion of proliferation, inflammation, and cancer.
The predicted metabolic profile is altered in patients with advanced neoplasia. A, Significant clustering by advanced neoplasia status on principal components analysis (PERMANOVA P = 0.001), with pronounced shifts in PC2 (P = 0.0003). B, Volcano plot demonstrating differentially abundant metabolites in advanced neoplasia. C, Significantly altered predicted levels of L-lactic acid, 2-ketobutyric acid, butyric acid, and L-tryptophan in advanced neoplasia. Plot capped at −4 for butyrate. P, Mann–Whitney test.
The predicted metabolic profile is altered in patients with advanced neoplasia. A, Significant clustering by advanced neoplasia status on principal components analysis (PERMANOVA P = 0.001), with pronounced shifts in PC2 (P = 0.0003). B, Volcano plot demonstrating differentially abundant metabolites in advanced neoplasia. C, Significantly altered predicted levels of L-lactic acid, 2-ketobutyric acid, butyric acid, and L-tryptophan in advanced neoplasia. Plot capped at −4 for butyrate. P, Mann–Whitney test.
Salivary microbiome data improves on clinical risk factors–based prediction of HGD and EAC
We next performed an exploratory analysis to determine whether salivary microbiome features in this cohort could be used to distinguish HGD/EAC from non-BE patients. As a baseline, we first trained a gradient-boosted decision trees model that uses clinical risk factors for EAC and tooth loss to classify HGD/EAC. The classifier was tested in cross-validation on patients not seen in the training of that model and achieved an AUROC of 0.84 (95% CI, 0.79–0.89). The same process was then used to train a classifier using microbiome data, whereas, within each training fold, 10 OTUs were selected based on a Kruskal–Wallis test. This classifier had an AUROC of 0.72 (95% CI, 0.65–0.79). Finally, a model trained on the combination of both microbiome data and EAC risk factors resulted in somewhat higher model accuracy, producing an AUROC of 0.88 (95% CI, 0.83–0.91; vs. clinical risk factor model AUROC 0.84, 95% CI, 0.79–0.89; DeLong P = 0.053; Supplementary Fig. S4).
Discussion
In this cross-sectional study of patients with and without BE, we detected marked shifts in the salivary microbiome with increased neoplasia, with changes that appeared to be most pronounced in patients with HGD or EAC. These changes included reduced diversity as well as significantly increased relative abundance of several taxa in the genus Streptococcus. Although tooth loss was more common in patients with HGD/EAC, we found that many of the salivary microbiome associations observed in BE and HGD/EAC persisted even when accounting for it.
Our findings add to the growing body of evidence that the oral microbiome is linked to the esophageal microbiome and is associated with esophageal neoplasia. In a case–control study of patients with EAC, BE, and controls, the EAC-associated microbiome had significantly reduced alpha diversity, similar to our observations in saliva (12). Interestingly, in that study 5/15 of the EAC tumors were dominated by Streptococcus spp. (relative abundance 69%–98%). In our salivary microbiome analyses, 4 of the 5 taxa most strongly associated with HGD/EAC were also Streptococcus spp. Our group conducted a randomized controlled trial and found that an antimicrobial mouth rinse can produce esophageal microbiome and tissue gene-expression changes, highlighting the relevance of oral bacteria to esophageal disease (14). In another cohort study analyzing mouth rinse samples from patients enrolled in two large cancer prevention studies, oral microbiome alterations were noted to precede an EAC diagnosis by several years (15). In a small study of 49 patients, we previously noted marked salivary microbiome alterations associated with BE and also with HGD/EAC (16).
Our study features the use of microbiome metabolic models to identify broad shifts in metabolites predicted to be produced by the saliva microbiome. Many of the predicted changes to metabolite outputs correspond with existing knowledge. Lactic acid, for example, was predicted to be increased in HGD/EAC. Lactic acid can serve as a major energy source for proliferative cancer cells and is known to activate hypoxia-inducible factors, which in turn contribute to proliferation, angiogenesis, and other neoplastic features (41). Our findings could therefore support the hypothesis that the oral and esophageal microbiota promotes EAC development and progression via production of metabolites (48). Our findings further correspond with a previous study detecting lactic-acid bacteria in many EACs (12). However, the biological significance of predicted metabolite production is unclear, and future studies are needed to validate these predictions and to elucidate the biological effects of specific bacterial metabolites on esophageal neoplasia.
Prior work has associated tooth loss and periodontal disease with increased risks of esophageal squamous cell cancer and gastric cancer (49, 50), and a recent study found an association between both tooth loss and periodontal disease and risk of EAC (18). Our results are in line with these recent analyses of data from the Nurses’ Health Study, which found an association between both tooth loss and periodontal disease and risk of EAC, but showed a decrease in association strength after adjusting for covariates (18). This indicates a potential confounding effect, as tooth loss is also associated with major alterations in oral microbiome composition (17). Our study offers an explanation for these associations, demonstrating that salivary microbiome composition is independently associated with HGD/EAC, even when adjusting for EAC risk factors and tooth loss. These findings suggest that the association between tooth loss and esophageal neoplasia is mediated by changes in the salivary microbiome, and that the salivary microbiome may represent a novel independent risk factor for EAC.
We performed exploratory analyses to assess whether the salivary microbiome could discriminate patients at the highest EAC risk. The salivary microbiome is highly suitable for diagnostics, as it is stable over time (51–53), especially compared with other body sites (54), and is resistant to perturbations (55). The addition of a microbiome-based classifier to EAC risk factors resulted in modest improvement in discrimination. However, the current study was not specifically designed to address this question, and future studies should explore further the salivary microbiome as a potential biomarker for HGD/EAC.
Important strengths of the current study include the relatively large sample size and the inclusion of oral health and hygiene information from patients. The large sample size allowed for the detection of significant microbiome alterations, even when correcting for multiple comparisons. Previous studies of the oral microbiome in BE and EAC have not included oral health and hygiene data, key potential confounders. The patients were well characterized, with data collected on key EAC risk factors including GERD history, BMI, and smoking, which permitted microbiome analyses adjusting for these variables. The BE patients in the study were demographically similar to BE populations from other studies, enhancing the generalizability of the findings. Lastly, novel methods for predicted microbiome metabolic profiling allowed for insights into functional correlates of the salivary microbiome alterations.
The study does have certain limitations. There were a relatively small number of nondysplastic BE patients, limiting analyses in this subgroup. A control group of patients undergoing a first upper endoscopy may have served as a better comparator than the current comparator group, many of whom had established upper gastrointestinal conditions. Analyses did not incorporate dietary intake; however, previous studies suggest that diet has minimal impact on salivary microbiome composition (56–58). Alcohol use data were not collected; whereas not an established risk factor for BE, it may have represented an important covariate to include in analyses. No conclusions can be drawn with regard to temporality in this cross-sectional study. It is possible that the observed salivary microbiome alterations were caused by BE-associated dysplasia or EAC, although we believe that this is unlikely. Many of the BE patients had previously undergone endoscopic eradication therapy and did not have residual BE. We found no differences in the salivary microbiome comparing those who had and had not been previously treated; however, future studies restricted to treatment-naïve patients are indicated.
Tooth loss was self-reported rather than measured, and periodontal disease was not directly assessed. Further, microbiome alterations associated with dental hardware could not be distinguished from tooth loss in this study. Community-scale metabolic models also have notable limitations. Our analysis was based on 16S rRNA gene sequencing, which does not allow us to tailor models to specific strains or genetic potential present in each sample. Additionally, although genome-scale models have been curated for common gut commensals, to our knowledge, such efforts have not been done for oral microbes. Consequently, some models may be missing, whereas existing ones may lack a representation of niche-specific metabolic capacity. Despite these limitations, these models allow a systematic application of biochemical and genetic knowledge to our analysis and raise interesting hypotheses that could be experimentally validated. The fungal mycobiome is increasingly being associated with gastrointestinal cancers but was not assessed as part of the current study.
In conclusion, patients with BE-associated HGD and EAC have a markedly altered salivary microbiome, and analyses of taxonomic alterations associated with stages from BE to EAC appear to indicate that these changes are most notable at the transition from LGD to HGD. Increased tooth loss was also observed in patients with HGD and EAC, although the salivary microbiome alterations were largely independent of tooth loss, suggesting that the association of tooth loss with HGD and EAC is mediated through the oral microbiome. There were marked increases in various taxa in the genus Streptococcus in HGD and EAC. In addition to the microbiome alterations, increased BE-associated neoplasia was associated with numerous changes to predicted bacterial metabolite production. Further work is warranted to identify the biological significance of the microbiome alterations, to establish causality, to validate metabolic shifts, and to determine whether they represent viable therapeutic targets for the prevention of progression in BE.
Authors' Disclosures
C.J. Lightdale reports grants from Erbe and Pentax and personal fees from Pinnacle outside the submitted work. J.A. Abrams reports other support from Exact Sciences, Castle Biosciences, and Pentax Medical outside the submitted work. No disclosures were reported by the other authors.
Authors' Contributions
Q.S. Solfisburg: Data curation, formal analysis, investigation, writing–original draft, writing–review and editing. F. Baldini: Formal analysis, investigation, writing–original draft, writing–review and editing. B. Baldwin-Hunter: Data curation, writing–review and editing. G.I. Austin: Formal analysis, writing–review and editing. H.H. Lee: Resources, data curation, writing–review and editing. H. Park: Data curation, formal analysis, writing–review and editing. D.E. Freedberg: Methodology, writing–review and editing. C.J. Lightdale: Resources, data curation, investigation, writing–review and editing. T. Korem: Conceptualization, resources, software, supervision, investigation, methodology, writing–review and editing. J.A. Abrams: Conceptualization, supervision, funding acquisition, methodology, writing–original draft, writing–review and editing.
Acknowledgments
This study was supported in part by the National Cancer Institute (U54 CA163004 and R01 CA238433 to J.A. Abrams) and the Digestive Disease Research Foundation (to Q.S. Solfisburg). T. Korem is a CIFAR Azrieli Global Scholar in the Humans and the Microbiome Program. D.E. Freedberg was supported in part by a Department of Defense Peer Reviewed Medical Research Program Clinical Trial Award (PR181960) and by a Columbia University Irving Scholar Award.
Note: Supplementary data for this article are available at Cancer Epidemiology, Biomarkers & Prevention Online (http://cebp.aacrjournals.org/).