Signatures within the esophageal microbiome are associated with host genetics, age, and disease

Abstract

Background: The esophageal microbiome has been proposed to be involved in a range of diseases including the esophageal adenocarcinoma cascade; however, little is currently known about its function and relationship to the host. Here, the esophageal microbiomes of 106 prospectively recruited patients were assessed using 16S rRNA and 18S rRNA amplicon sequencing as well as shotgun sequencing, and associations with age, gender, proton pump inhibitor use, host genetics, and disease were tested.

Results: The esophageal microbiome was found to cluster into functionally distinct community types (esotypes) defined by the relative abundances of Streptococcus and Prevotella. While age was found to be a significant factor driving microbiome composition, bacterial signatures and functions such as enrichment with Gram-negative oral-associated bacteria and microbial lactic acid production were associated with the early stages of the esophageal adenocarcinoma cascade. Non-bacterial microbes such as archaea, Candida spp., and bacteriophages were also identified in low abundance in the esophageal microbiome. Specific host SNPs in NOTCH2, STEAP2-AS1, and NREP were associated with the composition of the esophageal microbiome in our cohort.

Conclusions: This study provides the most comprehensive assessment of the esophageal microbiome to date and identifies novel signatures and host markers that can be investigated further in the context of esophageal adenocarcinoma development.

Keywords: Community types; Esophagus; Metagenomics; Microbiota; Single nucleotide polymorphisms.

Fig. 1

The esophageal microbiome clusters into different community types (esotypes). a Heatmap of relative abundances of the top 50 OTUs generated through 16S rRNA amplicon sequencing were used for a hierarchical cluster analysis (HCA) within MetaPhlan2. Dark blue to yellow correspond to 0.1–100% abundance. All available samples (n = 122) were utilized in this analysis. Taxonomy of OTUs is provided in Additional file 1.1. b Non-metric multidimensional scaling (nMDS) plot of Bray–Curtis resemblance generated from square root-transformed OTU relative abundances (all OTUs). OTU relative abundances were generated from 16S rRNA amplicon sequencing. Clusters from the HCA (panel a) were overlayed onto the nMDS plot. c PERMANOVA across HCA clusters of Bray–Curtis resemblance generated from square root-transformed OTU relative abundances. Test of the homogeneity of multivariate dispersions within groups at OTU level using PERMDISP showed no differences across clusters. ANOSIM generated a sample statistic of 0.42 and P = 0.001. d Relative abundances of the top 50 species generated through shotgun sequencing and MetaPhlan2 analysis were used for HCA. HCA was performed in MetaPhlan2. e nMDS plot of Bray–Curtis resemblance generated from square root-transformed species relative abundances (shotgun). All available shotgun samples were utilized in this analysis. Clusters from HCA of shotgun data (MetaPhlan2) were overlayed onto the nMDS plot. f PERMANOVA across shotgun HCA clusters of Bray–Curtis resemblance generated from square root-transformed species relative abundances. ANOSIM at the species level generated a sample statistic of 0.54 and P = 0.001

Fig. 2

Esophageal microbiome community types are defined by diversity and composition. a Principal component analysis of square root-transformed OTU relative abundances. The relative abundances of Haemophilus, Streptococcus, and Prevotella per subject were overlayed onto the PCA to define each cluster. Size of circle corresponds to relative abundance (%) of taxon. All available samples were utilized in this analysis. b Comparison analysis of phylum and genus relative abundances (%) generated from MEGAN6 according to the community types. Cluster 1, yellow; cluster 2, blue; cluster 3, red. Cluster 1 showed an enrichment of Prevotella and Haemophilus, cluster 2 showed an enrichment of Streptococcus, and cluster 3 showed an enrichment of Prevotella and Veillonella. c Correlations across species (shotgun MetaPhlan2) for each community type were calculated using SparCC and correlations greater than 0.2 or lower than − 0.2 were visualized using Cytoscape. Thickness of line reflects the strength of correlation and color reflects direction (green: positive; red: negative). A complete list of SparCC correlations within each cluster is provided in Additional file 1.9. d Alpha diversity measures for each community type. ANOVA with Tukey’s multiple comparison tests were used to calculate P values. Results related to species evenness is provided in Additional file 1.10

Fig. 3

Esophageal community types are functionally distinct. a nMDS plot of Bray–Curtis resemblance generated from square root-transformed KEGG pathway (level 3) relative abundances (generated using HUMAnN2). Clusters from HCA of shotgun data (MetaPhlan2) were overlayed onto the nMDS plot. All available samples were utilized in this analysis. b PERMANOVA across shotgun HCA clusters of Bray–Curtis resemblance generated from square root-transformed KEGG pathway relative abundances. ANOSIM at KEGG pathway level 3 generated a sample statistic of 0.46 and P = 0.001. c KEGG pathways identified using LEfSe analysis to be differentially abundant across each community type. All available samples within each cluster were utilized in this analysis. Blue, cluster 1; green, cluster 2; red, cluster 3. d MetaCyc pathways identified using LEfSe analysis to be differentially abundant across each community type. A full list of pathway names can be found in Additional file 1.14. All available samples within each cluster were utilized in this analysis. ANOSIM for MetaCyc pathways generated a sample statistic of 0.41 and P = 0.001. Blue, cluster 1; green, cluster 2; red, cluster 3

Fig. 4

Esophageal microbial signatures associated with the early stages of the esophageal adenocarcinoma cascade. a Microbial taxa identified using LEfSe analysis to be differentially abundant between GERD and subjects with a normal esophagus. The analysis was performed after stratifying the subjects according to community types. Green, normal; red: GERD. b Microbial taxa identified using LEfSe analysis to be differentially abundant between BE and subjects with a normal esophagus. Only samples designated as BE-Y (not BE-N or BE-GERD) were employed for this analysis. Red, BE. c Microbial taxa identified using LEfSe analysis to be differentially abundant between GM and subjects with a normal esophagus. Samples designated as GM-Y (not GM-N or GM-GERD) were employed for this analysis. Red, GM. d Correlations across species (shotgun MetaPhlan2) for each disease type were calculated using SparCC and correlations greater than 0.2 or lower than − 0.2 were visualized using Cytoscape. The thickness of line reflects the strength of the correlation, while color reflects direction (green, positive; red, negative). Samples designated as BE-Y were employed for this analysis. A complete list of SparCC correlations within each disease subgroup is provided in Additional file 3.8. e MetaCyc pathways identified using LEfSe analysis to be differentially abundant between disease type (GERD or BE) and subjects with a normal esophagus. Only samples designated as BE-Y were employed for this analysis. Green: normal; red: disease (GERD or BE)

Fig. 5

Presence of non-bacterial microbial taxa within the esophageal microbiome. a Relative abundance of viruses and fungi within the esophageal microbiome of each subject. Relative abundances were calculated using taxonomy arising from MEGAN6. This was used due to its capacity to detect microbial eukaryotes. b Relative abundances of specific eukaryotic and c viral taxa within each subject arising from the MEGAN6 analysis. Size of circle signifies the relative abundance levels of the organism. Size of circles ranges from 0.0094 (smallest) to 0.64% (largest) for viruses and 0.0097–1.08% for eukaryotes. Subjects are ordered by time of recruitment (left to right). The presence of Candida spp. and Saccharomyces were confirmed using 18S amplicon sequencing. We could not confirm the detection of Trichuris, Trichinella, and Loa loa; thus, these identifications should be taken with caution as they most probably arose from misclassifications within the metagenomic annotation process. No association with any of the clinical metadata was found

Fig. 6

Host genetic factors associated with the esophageal microbiome. a Host SNPs identified using MicrobiomeGWAS to be correlated with the Bray–Curtis resemblance matrix generated from square root-transformed species relative abundances (taxonomy arising from MetaPhlan2). Human SNPs were identified using the GATK toolkit onto the shotgun sequencing reads (depth of coverage (dp) = 2; minimum number of samples = 50). Blue line represents P = 0.1 and red line represents P = 0.05; all SNPs above the red line have significant P values. A complete list of SNPs across different thresholds is provided in Additional file 4.1–3. SNPs associated with microbiome composition mapped across most human chromosomes suggesting the analysis was not biased by low depth of coverage of the human genome. b PERMANOVA on Bray–Curtis resemblance matrix generated from square root-transformed species relative abundances (taxonomy arising from MetaPhlan2). Tests were applied across allele and genotype frequencies for human SNPs validated using Fluidigm custom SNPtype assays. A table of the genotyping results generated from the Fluidigm custom assays is provided in Additional file 4.4 and 5