Oral Microbiome Composition Reflects Prospective Risk for Esophageal Cancers

Esophageal cancer is the eighth most common cancer and sixth most common cause of cancer-related death worldwide (1). Because late-stage presentation is common in most cases, esophageal cancers are highly fatal; 5-year survival rates range from 15% to 25% in most countries (2). Consequently, there is a critical need for new avenues of prevention, risk stratification, and early detection.

The two main types, esophageal adenocarcinoma (EAC) and esophageal squamous cell carcinoma (ESCC), differ greatly in incidence, geography, and etiology. ESCC, the most common type worldwide, predominates in developing countries, while EAC has become the predominant type in developed countries as incidence rates continue to rise (2, 3). Known risk factors include gastro-esophageal reflux disease (GERD), obesity, low fruit/vegetable intake, and smoking for EAC, and alcohol drinking, low fruit/vegetable intake, and smoking for ESCC (4), but the etiology of these diseases cannot be fully explained by these factors.

Recently, upper digestive tract microbiota have been suggested to play a role in esophageal cancer etiology, and in particular in the rising incidence of EAC in developed countries (5). The complex microbial community of the upper digestive tract, consisting of mutualists, commensals, and pathogens, could facilitate carcinogenesis via activation of Toll-like receptors (6), or protect against carcinogenesis via synthesis of vitamins or providing barriers to pathogen invasion (5). Cross-sectional studies report distinct differences in upper digestive tract microbiota between GERD (7–9), Barrett's esophagus (an EAC precursor; refs. 7–10), EAC (7, 11), esophageal squamous dysplasia (ESD, an ESCC precursor; ref. 12), or ESCC (13) cases and controls. In addition, periodontitis (a disease of oral dysbiosis) may be associated with increased esophageal cancer risk (14). However, no studies have prospectively examined whether upper digestive tract microbiota influence risk for subsequent esophageal cancer.

We hypothesized that oral microbiota influence development of esophageal cancer. The oral microbiota shape the esophageal microbiome (15), due to migration of oral bacteria to the esophagus (16) and, therefore, may contribute to esophageal carcinogenesis. We conducted a prospective study nested in two large U.S. cohorts, to determine whether oral microbiota are associated with subsequent EAC or ESCC risk.

Participants were drawn from two U.S. cohorts: the NCI Prostate, Lung, Colorectal, and Ovarian (PLCO) Cancer Screening Trial cohort and American Cancer Society (ACS) Cancer Prevention Study II (CPS-II) Nutrition cohort. Characteristics of these cohorts are comparable, with both cohorts collecting oral wash samples and comprehensive demographic information, and following prospectively for cancer incidence.

PLCO (17) is a large population-based randomized trial designed to determine effects of screening on cancer-related mortality in men and women aged 55 to 74, recruited in 1993 to 2001, and followed for cancer incidence. Participants were randomized to a screening or control arm. Oral wash samples were collected in the control arm only (n = 52,000). Incident cancers were ascertained by annual mailed questionnaire and verified through medical records or death certificates.

CPS-II (18) includes >184,000 participants, aged 50 to 74 who completed a mailed baseline questionnaire in 1992. Follow-up questionnaires have been sent to cohort members every other year to update information and ascertain incident cancers, which are also verified through medical records, state registries, or death certificates. During 2000 to 2002, oral wash samples were collected from 70,004 participants.

Incident cases were cohort participants diagnosed with esophageal cancer any time after oral wash collection (collection to diagnosis time ranged from <1 year to 9 years; first quartile, median, third quartile = 1, 3, 5 years) and had no prior cancer history (except nonmelanoma skin cancer). Matched controls were selected at a case:control ratio of 1:2 by incidence density sampling without replacement among participants who provided an oral wash sample in the same year as the index case, had no cancer at or prior to index case diagnosis, and were of the same cohort, age, sex, and race as the index case.

A total of 368 samples were provided and successfully sequenced, including 117 complete sets (1 case: 2 controls), 2 reduced sets (1 case: 1 control), and 13 unmatched controls (due to missing case, case failing sequencing, or nonesophageal case). On the basis of ICD morphology codes (EAC: 8140, 8144, 8480, 8481, 8560; ESCC: 8070, 8071, 8072, 8074, 8052), we included 81 EAC cases (with 160 matched controls) and 25 ESCC cases (with 50 matched controls) in the current analysis (N = 316). Cases of other or missing morphology (n = 13), their matched controls (n = 26), and unmatched controls (n = 13) were excluded.

This study was conducted in accordance with the U.S. Common Rule and approved by the IRB of New York University School of Medicine (New York, NY), NCI (Bethesda, MD), and ACS, and participants provided informed consent.

Covariate information was extracted from questionnaires preceding oral sample collection for each participant. Body mass index (BMI) was categorized as normal or underweight (BMI < 25 kg/m²), overweight (25 ≤ BMI < 30 kg/m²), or obese (BMI ≥ 30 kg/m²). Smoking status was classified as never, former, or current. Drinking level was classified as never, moderate, or heavy (19). Servings of fruits and vegetables per day, derived from food frequency questionnaire responses, were categorized as low or high based on cohort-specific medians.

Participants were asked to swish with 10 mL Scope mouthwash (P&G) and expectorate into a tube (17, 18). Samples were shipped to each cohort's biorepository and stored at −80°C. The oral microbiome is highly stable over time (20–22) and shows much greater interindividual than intraindividual variation, indicating that a one-time oral sample collection is appropriate for assessing oral microbial risk factors in a cohort study.

We extracted DNA from oral wash samples using the PowerSoil DNA Isolation Kit (Mo Bio). Barcoded amplicons were generated covering the 16S rRNA gene V4 region using F515/R806 primers. The PCR reaction used FastStart High Fidelity PCR System, dNTP pack (Roche) as follows: initial denaturing at 94°C for 3 minutes, followed by 25 cycles of 94°C for 15 seconds, 52°C for 45 seconds and 72°C for 1 minute, and a final extension at 72°C for 8 minutes. PCR products were purified using Agencourt AMPure XP (Beckman Coulter Life Sciences), quantified using Agilent 4200 TapeStation (Agilent Technologies), pooled at equimolar concentrations and sequenced on Illumina MiSeq with a 300-cycle (2 × 151 bp) reagent kit.

Paired-end reads were joined and demultiplexed, and poor-quality reads excluded, using default parameters in QIIME (23). The 11,422,831 quality-filtered reads (from N = 368 samples) were clustered into operational taxonomic units (OTU) against the Human Oral Microbiome Database (HOMD) reference sequence collection (version 14.5; ref. 24), and assigned HOMD taxonomy, using QIIME script pick_closed_reference_otus.py (23). This method discards reads not matching the database, leaving 11,074,719 reads [mean ± SD = 30,094 ± 21,059; range = (4,965–203,242)] and 569 OTUs. We generated a phylogenetic tree from aligned HOMD reference sequences using FastTree (25).

All samples underwent DNA extraction and sequencing in the same laboratory, with personnel blinded to case/control status. DNA from volunteer oral wash samples was included in the sequencing batches: six replicates from each of 4 volunteers in the CPS-II batch, and eight replicates from each of the same 4 volunteers in the PLCO batch. Intraclass correlation coefficients for the Shannon diversity index and relative abundance of major oral phyla were high (Supplementary Table S1), and principal coordinate analysis of UniFrac distances (26) showed clustering of repeat samples for each volunteer, indicating excellent reproducibility (Supplementary Fig. S1).

We used multiple imputation (“mice” package, R; ref. 27) to impute missing data for three important predictors of esophageal cancer, BMI, alcohol drinking, and fruit and vegetable intake. A total of 23 participants (7.3%) were missing BMI, 36 (11.4%) were missing alcohol drinking, and 39 (12.3%) were missing fruit and vegetable intake (% missing by case/control group shown in Table 1). Predictors of BMI category (<25, 25–30, ≥30) used in imputation were sex, race, age, cohort, smoking status, education level, and ethanol intake. Predictors of alcohol drinking (none, moderate, heavy) and fruit/vegetable intake (low or high) used in imputation were sex, race, age, cohort, smoking status, education level, and continuous BMI. Ten imputed datasets were used in analysis, and we present pooled estimates and P values.

α-Diversity (within-subject diversity) was assessed by richness and the Shannon diversity index, calculated in 100 iterations of rarefied OTU tables of 4,500 sequence reads per sample. This depth was chosen to sufficiently reflect sample diversity (Supplementary Fig. S2) while retaining all participants. We examined whether α-diversity differed between cases and controls in conditional logistic regression using matched sets as strata and adjusting for smoking status, BMI category, and alcohol drinking level.

β-Diversity (between-subject diversity) was assessed at OTU level using unweighted and weighted UniFrac distances (26). Permutational multivariate analysis of variance (PERMANOVA; “adonis” function, “vegan” package, R; ref. 28) was used to examine statistically whether overall bacterial community composition differed by case/control status, using matched sets as strata and adjusting for smoking status, BMI category, and alcohol drinking level.

The 569 OTUs were agglomerated to 12 phyla, 26 classes, 42 orders, 70 families, 149 genera, and 513 species. We applied the centered log-ratio (clr) transformation (29) to the taxa counts at each level (e.g., phylum, class, etc.) after adding a pseudocount of 1. We used conditional logistic regression, using matched sets as strata and adjusting for smoking status, BMI category, and alcohol drinking level, to determine whether abundance of bacterial taxa predicts esophageal cancer risk. This analysis included only taxa present in ≥15% of the 316 participants (10 phyla, 20 classes, 28 orders, 46 families, 85 genera, 266 species), to exclude rare taxa and thereby minimize the number of statistical tests conducted. A priori species of interest were “red complex” periodontal pathogens: Tannerella forsythia (T. forsythia), Porphyromonas gingivalis (P. gingivalis), and Treponema denticola (T. denticola; ref. 30). For other taxa, P values were adjusted for the FDR.

Metagenome content was predicted using PiCRUSt (31). Because PiCRUST gene content is precomputed for the GreenGenes database of 16S rRNA genes, for this analysis, we performed closed-reference OTU picking against the GreenGenes database prior to PiCRUST. The 5507 KEGG (32) gene orthologs were grouped into 270 KEGG pathways. We applied the clr transformation (29) to pathway counts after adding a pseudocount of 1, filtered to include pathways present in ≥15% of participants (255 pathways), and used conditional logistic regression, as described above, to determine whether abundance of functional pathways predicts esophageal cancer risk.

Ecological networks among species were inferred using the SPIEC-EASI (SParse InversE Covariance Estimation for Ecological Association Inference) algorithm (33). This statistical method, designed for ecological network inference from amplicon sequencing datasets, accounts for compositional data structure using the clr transformation and assumes a sparse underlying ecological association network. We applied SPIEC-EASI separately to EAC cases and matched controls, and ESCC cases and matched controls. The “igraph” package in R was used for network visualization.

All statistical tests were two-sided. A P < 0.05 was considered of nominal significance, and an FDR-adjusted P value (q-value) <0.10 was considered significant after multiple comparisons adjustment. Analyses were conducted using R 3.2.1.

Cases and their matched controls did not differ on matching factors (Table 1). Although obesity, low fruit/vegetable intake, and smoking are recognized risk factors for EAC, and alcohol drinking, low fruit/vegetable intake, and smoking are recognized risk factors for ESCC, only alcohol drinking was associated with ESCC (P = 0.004).

EAC and ESCC cases did not differ significantly from matched controls in oral α-diversity, as measured by species richness and the Shannon diversity index, or overall oral microbiome composition (β-diversity), as measured by unweighted and weighted UniFrac distances (Supplementary Table S2).

For the a priori “red complex” periodontal pathogens (30), a doubling of T. forsythia abundance relative to the geometric mean of all taxa was associated with 1.21 [95% confidence interval (CI), 1.01–1.46] times higher odds of EAC (P = 0.04), while abundance of P. gingivalis and T. denticola was not associated with EAC risk (Table 2; Fig. 1). We identified several other oral taxa nominally associated with EAC risk (Table 3; Fig. 1), although none reached the significance threshold after FDR adjustment (all q-value > 0.32). Increased abundance of species Actinomyces cardiffensis, Selenomonas oral taxon 134, and Veillonella oral taxon 917 was associated with higher EAC risk (all P < 0.05). Conversely, increased abundance of Corynebacterium durum, Prevotella nanceiensis, Streptococcus pneumoniae, Lachnoanaerobaculum umeaense, Oribacterium parvum, Solobacterium moorei, Neisseria sicca, Neisseria flavescens, and Haemophilus oral taxon 908 was associated with lower EAC risk (all P < 0.05). Additional adjustment for fruit/vegetable intake did not impact effect estimates (percent change in β-coefficient for all nominally significant taxa <12%).

We observed that the majority of these species were associated with each other in an ecological network analysis (Fig. 2A). The protective species in phylum Proteobacteria (Neisseria sicca, Neisseria flavescens, and Haemophilus oral taxon 908) were closely connected, as were some of the protective species in phylum Firmicutes (Solobacterium moorei, Oribacterium parvum, and Lachnoanaerobaculum umeaense). Some of the species formed their own networks (i.e., unrelated to other EAC-associated species), including Streptococcus pneumoniae and Selenomonas oral taxon 134.

We additionally explored heterogeneity of taxon abundance–EAC associations by years from oral wash collection to diagnosis (≤ or > median of 3 years), cohort (CPS-II or PLCO), smoking status (ever or never), obesity (nonobese or obese), and fruit/vegetable intake (low or high). Taxon findings were consistent across years to diagnosis subgroups (all P_interaction > 0.12; Supplementary Table S3). Similarly, taxon findings were largely consistent across cohorts (Supplementary Table S4; Supplementary Fig. S3); in particular, Streptococcus pneumoniae, Solobacterium moorei, Veillonella oral taxon 917, Neisseria sicca, Neisseria flavescens, and Haemophilus oral taxon 908 showed homogenous associations with EAC in both cohorts. Selenomonas oral taxon 134 was associated with higher EAC risk in the PLCO cohort only (P_interaction = 0.02). Similarly, the periodontal pathogens tended to be associated with higher EAC risk only in PLCO (P_interaction = 0.11, 0.35, and 0.04 for P. gingivalis, T. forsythia, and T. denticola, respectively). When stratifying by smoking status, we observed that Lachnoanaerobaculum umeaense was associated with lower EAC risk only in smokers (P_interaction = 0.02; Supplementary Table S5), while other taxon–EAC associations did not differ significantly between ever and never smokers (P_interaction > 0.19). When stratifying by obesity, we observed that Actinomyces cardiffensis was associated with higher EAC risk only in nonobese (P_interaction = 0.02), while other taxon–EAC associations did not differ significantly between nonobese and obese (P_interaction > 0.11; Supplementary Table S6). Finally, when we stratified by fruit and vegetable intake, order Actinomycetales was associated with higher EAC risk only in those with higher fruit and vegetable intake (P_interaction = 0.05), while other interactions were nonsignificant (P_interaction > 0.18; Supplementary Table S7).

The periodontal pathogen P. gingivalis was marginally associated with higher ESCC risk [OR (95% CI) = 1.30 (0.96–1.77); P = 0.09; Table 2; Fig. 1]. Several other species were nominally associated with ESCC risk (Table 3; Fig. 1), although none reached the significance threshold after FDR adjustment (all q-value > 0.80). Increased abundance of Prevotella nanceiensis, Bergeyella oral taxon 322, Neisseria weaveri, and Treponema vincentii was associated with higher ESCC risk, while increased abundance of Prevotella oral taxon 306 and Aggregatibacter paraphrophilus was associated with lower ESCC risk (all P < 0.05). Additional adjustment for fruit/vegetable intake did not impact effect estimates (percent change in β-coefficient for all nominally significant taxa <11%). We did not perform stratified analysis of taxonomic findings for ESCC due to small sample size.

All of the species nominally associated with ESCC were associated with each other in an ecological network analysis (Fig. 2B). Interestingly, Treponema vincentii, which was associated with increased ESCC risk and has been previously associated with periodontal disease (34, 35), was linked to other periodontal pathogens (P. gingivalis, T. forsythia) in the ecological network.

Analysis of inferred metagenomes revealed a number of metabolic pathways nominally associated with EAC risk (Table 4), although none reached the significance threshold after FDR adjustment. Increased abundance of endocytosis, sulfur relay system, biosynthesis of siderophore groups, and bisphenol degradation pathways was associated with higher EAC risk, and α-linolenic acid (ALA) metabolism and carotenoid biosynthesis pathways with lower risk (all P < 0.05). We did not identify any pathways associated with ESCC risk. Species Neisseria sicca and Neisseria flavescens, associated with reduced EAC risk, were positively correlated with the protective carotenoid biosynthesis and ALA metabolism pathways (Fig. 3).

In this first prospective study of oral microbiota and esophageal cancer risk, we did not observe significant associations between overall microbiota diversity or composition and subsequent EAC or ESCC risk. However, several species were nominally associated with risk, among them potential pathogens and also commensal species. Interestingly, bacterial taxon associations observed were unique to either EAC or ESCC, in line with the fundamentally different origins of these cancer types. We also show replication of several taxonomic findings in both the CPS-II and PLCO cohorts. Our biologically plausible findings warrant further investigation in larger studies, to fully explore prospects of modulating the oral microbiota for esophageal cancer prevention or utilizing it for risk stratification and early detection.

Studies of oral disease and cancer provide evidence that oral health (tooth loss, poor oral hygiene, and possibly periodontal disease) is linked to esophageal cancer risk (14, 36–38). We observed that T. forsythia was associated with higher EAC risk, and P. gingivalis with ESCC risk. These two species are members of the “red complex” of periodontal pathogens, that is, the species most strongly associated with severe periodontitis (30). A recent report revealed that P. gingivalis was detected at a higher rate in ESCC tumor tissue, compared with adjacent normal and healthy control mucosa; moreover, P. gingivalis presence was associated with ESCC lymph node metastasis and decreased survival time (39). More research is needed to determine whether periodontal disease and/or periodontal pathogens play a role in EAC/ESCC carcinogenesis, particularly as periodontal pathogen–EAC risk associations were inconsistent between the CPS-II and PLCO cohorts.

Several small studies have characterized the esophageal microbiota in relation to EAC (7, 11) or its precursors, GERD (7–9) and Barrett's esophagus (7–10). Campylobacter species were shown to dominate GERD and Barrett's esophagus biopsies compared with controls in two culture-based studies of subjects from the United Kingdom (7, 10). Yang and colleagues surveyed 16S rRNA genes from distal esophageal biopsies of 12 controls, 12 GERD patients, and 12 Barrett's esophagus patients in the United States (8); they observed a distinctly different microbial composition in GERD and Barrett's esophagus patients compared with controls, characterized by greater diversity, decreased Streptococcus, and increased abundance of Gram-negative anaerobes, including Veillonella, Neisseria, Prevotella, Campylobacter, Porphyromonas, Fusobacterium, and Actinomyces. Similarly, Japanese patients with Barrett's esophagus had decreased Streptococcus and increased Veillonella, Neisseria, and Fusobacterium in distal esophageal biopsies compared with controls (9). Finally, Zaidi and colleagues observed decreased abundance of Streptococcus pneumoniae in dysplastic, tumor-adjacent normal, and EAC biopsy samples compared with normal and Barrett's esophagus samples from U.S. patients (11). We observed an inverse association between Streptococcus pneumoniae and incident EAC, consistent with above-mentioned studies. In contrast to above-mentioned studies, we observed an inverse association of genus Neisseria with EAC risk. Neisseria species are oral cavity commensals (40), and we and others previously showed that oral Neisseria are depleted by cigarette smoking (41–43), a cause of EAC. Interestingly, we found that Neisseria were only associated with lower EAC risk in smokers (although interaction was not significant), possibly suggesting a joint effect of smoking and Neisseria depletion. Differing findings from previous literature may relate to differences in study design (cross-sectional vs. prospective) and sample origin (biopsy vs. oral).

Other studies have characterized the microbiota related to ESCC (13) and its precursor, ESD (12). Yu and colleagues observed that lower microbial richness and altered composition of upper digestive tract microbiota were associated with ESD in Chinese subjects (12). Likewise, Chen and colleagues reported differences in carriage and/or relative abundance of oral genera between 87 ESCC cases and 85 controls, including increased relative abundance of Prevotella, Streptococcus, and Porphyromonas in ESCC cases. These authors did not report findings at species level, making comparison with our mostly species-level findings for ESCC difficult.

Analysis of inferred metagenomes revealed several pathways associated with EAC, albeit not after FDR adjustment; some appeared biologically plausible. Bacterial carotenoid biosynthesis was associated with lower EAC risk, with Neisseria species potentially contributing to this protective pathway. Carotenoids are phytochemicals in fruits and vegetables, many acting as antioxidants (44). Higher fruit and vegetable intake and higher β-carotene intake have been associated with reduced EAC risk (4, 45). In addition, β-carotene therapy was shown to ameliorate GERD symptoms (46). Bacterial biosynthesis of siderophores (iron-chelating compounds) was associated with higher EAC risk. Although excessive iron may promote carcinogenesis (47) and iron chelation has been considered as a potential EAC therapy (48), iron is an essential trace element with deficiency leading to inflammation (49). Bacterial siderophore synthesis may upset iron homeostasis and thus might increase EAC risk. These inferred metagenomic functions provide insight into bacterial actions that may potentially impact EAC risk and warrant further investigation with full metagenomic sequencing.

Strengths of our study included the prospective design, comprehensive 16S rRNA gene sequencing, inclusion of two cohorts, and adjustment for EAC/ESCC risk factors throughout analysis. Our study also had several limitations. Lack of periodontal status of participants did not allow us to determine whether periodontal pathogens are implicated independently of periodontal disease. We also lacked data on presence of esophageal cancer precursor conditions (i.e., GERD and Barrett's esophagus) in the participants, which could mediate or confound oral microbiome-esophageal cancer associations and data on medications (e.g., proton-pump inhibitors, antibiotics), which could confound these associations (50, 51). In addition, although our study is the largest of its kind, case sample sizes (n = 81 EAC and n = 25 ESCC) remained small, limiting statistical power to detect FDR-adjusted significant associations, and our study population was mostly white, limiting generalizability.

In summary, we found evidence that specific bacterial pathogens may play a role in esophageal cancer risk, whereas other bacterial types may be associated with reduced risk. Larger studies are needed to confirm our findings, particularly among smokers and nonsmokers to clarify joint effects, followed by experimental animal models to clarify causal relationships. Identification of oral bacteria causal or protective in esophageal cancer could lead to interventions for their eradication or colonization in at-risk individuals. Continued study of oral microbiota in esophageal cancer may lead to actionable means for prevention of this highly fatal disease.

No potential conflicts of interest were disclosed.

Z. Pei is a staff physician at the Department of Veteran Affairs New York Harbor Healthcare System. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH, the U.S. Department of Veterans Affairs, or the United States Government.

Conception and design: Z. Pei, R.B. Hayes, J. Ahn

Development of methodology: J. Wu, Z. Pei, L. Yang, M.P. Purdue, R.B. Hayes, J. Ahn

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.P. Purdue, E.J. Jacobs, S.M. Gapstur, R.B. Hayes, J. Ahn

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): B.A. Peters, N.D. Freedman, E.J. Jacobs, S.M. Gapstur, R.B. Hayes, J. Ahn

Writing, review, and/or revision of the manuscript: B.A. Peters, Z. Pei, M.P. Purdue, N.D. Freedman, E.J. Jacobs, S.M. Gapstur, R.B. Hayes, J. Ahn

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): Z. Pei, N.D. Freedman, J. Ahn

Study supervision: Z. Pei, J. Ahn

This work was supported in part by the NCI at the NIH (R01CA159036, U01CA182370, R01CA164964, R03CA159414, P30CA016087, and R21CA183887). The NYU School of Medicine Genome Technology Center where samples were sequenced is partially supported by the Perlmutter Cancer Center support grant P30CA016087. The American Cancer Society (ACS) funds the creation, maintenance, and updating of the Cancer Prevention Study II cohort.

Samples were sequenced at the NYU School of Medicine Genome Technology Center.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.