Gastroesophageal reflux disease (GERD) is one of the most common gastrointestinal tract (GIT) disorders, with prevalence reaching up to 20% in Western countries.^1,2 GERD complications include an inflammation of the esophagus, known as reflux esophagitis (RE), Barrett´s esophagus (BE), and esophageal adenocarcinoma (EAC);³ however, only a minority of patients with GERD develop BE or EAC. The long-term exposure of the normal squamous epithelium to gastroduodenal fluid results in an inflammation of the esophagus and, subsequently, in replacement of squamous epithelium with columnar epithelium (intestinal metaplasia) in the lower esophagus.⁴ Intestinal metaplasia can then lead to the development of BE, which is acknowledged as a pre-cancerous lesion for EAC development; the exact mechanism and conditions of the conversion in this inflammation-related carcinogenesis have, however, not been elucidated yet.⁵

The pathogenesis of GERD includes multiple external and internal exposure factors, including obesity, lack of exercise, unhealthy lifestyle (alcohol abuse, smoking, etc.),⁶ and genetic predispositions.^1,7 Most patients with GERD, especially those with EAC, are older men.⁸ Besides these factors, esophageal microbiota has also been associated with GERD pathogenesis.^9,10 While the quantity of microbiota in the healthy esophagus is lower than in healthy oral mucosa, the composition of these sites is very similar,¹¹ with the predominance of grampositive aerobic bacterial genera, especially by the genus Streptococcus. Specimens from the esophagus with GERD-induced pathological changes have, however, been associated with higher relative abundances of gramnegative anaerobes or microaerophiles, such as Veillonella, Prevotella, Haemophilus, Neisseria, Rothia, Granulicatella, Campylobacter, Porphyromonas, Fusobacterium, and Actinomyces.^12–14 In this context, the role of pharmaceuticals, such as proton pump inhibitors (PPIs) that are commonly used by these patients for reducing stomach acid and/or Helicobacter pylori eradication (in combination with antibiotics), must also be considered, as they can influence the GIT microbiome diversity and composition.^15,16

Although several bacteria were associated with GERD and its complications, scientific literature still lacks a comprehensive description of the GIT bacteriome (i.e., complete genetic information of bacteria in the entire GIT), which could shed light on the interconnectedness of different parts of the GIT in these patients. Identification of bacteriome features specific to these patient groups could help us to better understand the pathogenesis of the disease and the relationships among GERD complications. Here, we aimed to describe, compare, and correlate bacteriomes in biopsies from macroscopically pathological esophageal tissues and those from other parts of the GIT (oral and rectal swabs, gastric body and gastric antrum biopsies, and duodenal biopsies) in patients with RE, BE, and EAC. In addition, as several studies have suggested the possible use of some oral bacteriome characteristics as non-invasive predictive biomarkers for BE/EAC,^17–19 we also aimed to address this hypothesis by comparing relative abundances of bacterial genera in oral and rectal swabs taken from patients with RE/BE/EAC.

One part of the biopsy from each site was processed to formalin-fixed paraffin-embedded (FFPE) blocks and examined by pathologists for histological confirmation of RE, BE, or EAC (patients with negative or inconclusive results were excluded from the study); standard immunohistochemical analysis for H. pylori was performed. The other part of the biopsy from each site, as well as the oral and rectal swabs, were placed in RLT buffer (from the AllPrep DNA/RNA Mini Kit; QIAGEN, Germany) and immediately stored at -80°C until the extraction of nucleic acids.

DNA was isolated from all samples (including 46 negative controls, NCs, and internal standard) using the AllPrep DNA/RNA Mini Kit (QIAGEN, Germany). In all matrices, the total content of bacterial DNA was established using qPCR targeting the 16S rRNA region. Subsequently, all samples and NCs were spiked with an internal standard and analyzed by 16S RNA amplicon sequencing on the Illumina platform. Details on laboratory and bioinformatic metagenomic analyses were described in our previous work.²⁷

Statistical analyses were carried out in R (v4.1.2). A multi-step filtering of samples was performed prior to the analysis. A threshold of 5,000 sequencing reads was established to prevent analysis of samples with an insufficient sequencing depth. Next, all reads originating from the internal standard were removed from the dataset, and every sample was tested for dissimilarity in bacteriome composition against NCs. Additional exclusion criteria based on the quality of the obtained data are shown in the study workflow in Fig. 1.

The genera present in ≥ 20 samples from any single GIT location (or NCs) were analyzed using the ALDEx2 framework. Clustering analysis of bacterial genera was also performed to identify the similarity between individual GIT sites.

Sequencing depth and alpha diversity, including the number of distinct amplicon sequence variants (ASVs) and Shannon diversity index, were then calculated on the filtered samples using the vegan package. For all downstream analyses, counts were collapsed to the genus level; any ASVs that could not be assigned a genus were grouped under the category “other”. An overview table of clinical variables was generated using an on-site tool based on the gt package.

The ratios of relative abundances of anaerobic to aerobic bacterial genera (anaerobes/aerobes ratio) and gramnegative to grampositive bacterial genera (G^-/G⁺ ratio) were computed based on all determined genera for which both statuses (i.e., aerobicity and Gram staining) were known, using, on average, 82% of the total bacterial DNA present in the sample. Dissimilarity in bacteriome composition was calculated as the sum of differences in relative abundances weighted by the difference in quantiles across all genera. This dissimilarity index was used in non-metric multidimensional scaling (nMDS) and to test differences in composition between non-NC samples and negative controls.

Comparisons of continuous variables between study groups, such as the total content of bacterial DNA and alpha diversity metrics, were carried out using two-sample Wilcoxon rank-sum tests, performed separately for each sample matrix. P-values from these tests were adjusted for multiple comparisons within the sample matrix by the Benjamini–Hochberg false discovery rate procedure. Results were visualized as boxplots in R using ggplot2, with styled significance annotations added via ggtext. Adjusted p-values are indicated by asterisks: * for P < .05, ** for P < .01, and *** for P < .001.

We used ALDEx2 to assess bacteriome composition differences. Read counts were transformed using the central log-ratio transformation using ALDEx2::aldex.clr with 256 Monte Carlo replicates and the scale parameter gamma = 0.25. To detect genera present in the analyzed matrices, a general linear model was built with NCs serving as a baseline. Only genera with significantly higher CLR-transformed relative abundances in the given matrix compared to NCs were considered as present and used in further analysis. Next, a series of Wilcoxon rank-sum tests was performed for every sample matrix, comparing study groups using the ALDEx2::ttest. The beta diversity was visualized using t-distributed stochastic neighbor embedding (t-SNE) with the Rtsne library.

Pairwise Pearson correlations were calculated on CLR-transformed abundance data, and p‐values were adjusted for multiple testing using the Benjamini–Hochberg procedure. Correlations with an absolute coefficient greater than 0.4 and an adjusted P‐value below .05 were retained and visualized as a network in R using the igraph package.

The total content of bacterial DNA in oral samples of the BE group was significantly higher than in the RE (P_ADJ < .05) and EAC (P_ADJ < .001, Fig. 2A) groups. In esophageal samples, the total content of bacterial DNA increased with the severity of GERD complications (P_ADJ < .001 for both RE-EAC and BE-EAC contrasts). It was also higher in duodenal samples from patients with EAC than with BE. The analysis of sequencing depth (Fig. 2B) revealed a similar pattern in esophageal samples – the sequencing depth was increased in EAC samples compared to groups with other GERD complications (P_ADJ < .01 for RE-EAC and < .05 for BE-EAC). The same result was observed in the gastric antrum biopsy (P_ADJ < .05 for both RE-EAC and BE-EAC). The esophagus was the only matrix in which a significant difference in bacteriome richness (measured as the number of distinct ASVs) was detected between GERD complications, with EAC samples containing more distinct ASVs than RE samples (P_ADJ < .01, Fig. 2C). An even stronger result was obtained in Shannon diversity (bacteriome evenness), where esophageal samples from patients with EAC showed higher diversity than both RE (P_ADJ < .001) and BE (P_ADJ < 0.05) samples, see Fig. 2D. The anaerobes/aerobes ratio (Fig. 2E) was significantly higher in oral swabs from patients with BE compared to both RE (P_ADJ < .05) and EAC (P_ADJ < 0.01) patient groups. Finally, the G^-/G⁺ ratio was significantly higher in esophageal samples from patients with BE compared to those with RE (P_ADJ < 0.05). In gastric body biopsies, this ratio was significantly higher in samples from patients with BE than with EAC (p_ADJ < 0.05), see Fig. 2F.

The overall bacterial composition is visualized using t-SNE in Fig. 3, showing a distinct cluster of rectal swabs and separation of oral samples from the rest of the GIT. The esophageal, gastric, and duodenal biopsy samples were overlapping and, in two dimensions, unseparable. Supplementary Table S1 shows the ten bacterial genera (based on their median relative abundance) that were most abundant in each of the six GIT sites. Genera Streptococcus, Veillonella, Haemophilus, and Prevotella are among the five most abundant genera in all GIT samples except rectal swabs (although Prevotella is among the top five genera even in the latter matrix).

Figure 4 shows an MDS plot of mean dissimilarity in bacteriome composition between individual GIT sites based on pairwise testing of esophageal samples with different severities of GERD complications and individual-matched samples from other locations. Note that in this plot, the mutual distances of esophageal samples for individual diagnoses do not indicate their true (dis)similarity as their mutual pairwise testing was, obviously, impossible. In esophageal tissues, bacteriome compositions of the EAC and RE were relatively similar to oral samples, while high dissimilarity was observed between the BE and oral samples. Bacteriome composition of duodenal samples from all three study groups was relatively similar to that of esophageal samples. Rectal swabs were omitted from this analysis due to the very low intersection of bacterial genera with other GIT sites (leading to a high dissimilarity to all other tissues, preventing a meaningful depiction).

The ALDEx2 framework was applied independently to each sample matrix. Across all matrices, a total of 27 bacterial genera/sample matrix combinations yielded at least one statistically significant contrast. All selected genera and their CLR-transformed abundance profiles across matrices and study groups are presented in Supplementary Fig. 2.

The highest number of significant contrasts between study groups for individual genera (12 genera, 18 contrasts in total) was found in oral samples. Porphyromonas was the only genus showing significant contrasts for both RE-EAC and BE-EAC in oral swabs; the same can be said of Veillonella in gastric body as well as gastric antrum samples. Significant contrasts were found also for both RE-EAC and RE-BE comparisons for Streptococcus in esophageal samples. In rectal swabs, differences between patients with EAC and RE were found only for Faecalibacterium. For full results including possible species classification, see Table 2 and Supplementary Figure S3, depicting distributions and median values of CLR-transformed relative abundances of genera with at least one significant contrast (p < .05).

The association of several bacteriome characteristics in multiple parts of the GIT with different severity of GERD complications may be caused by an interplay of the pathological mucosal changes and the bidirectional distribution of bacteria through the entire GIT (both due to swallowing and gastroesophageal reflux).¹⁴ In all, 923 samples collected from the entire GIT of 158 patients were analyzed in this study, providing a comprehensive picture of the bacteriomes in each patient, while also enabling us to compare findings from each location. We found that some bacteriome characteristics in biopsies from macroscopically pathological esophageal tissues differed among groups of patients with RE, BE, and EAC. Moreover, strong correlations were detected for some bacterial genera, such as Prevotella, throughout the entire GIT. More interestingly, high or low relative abundances of some bacterial genera, such as Porphyromonas on the oral mucosa and Faecalibacterium in fecal samples, were associated with EAC.

Our study identified similarities and differences among multiple GIT sites and groups of patients with RE, BE, and EAC. A gradual trend in several parameters (increase in the total content of bacterial DNA and alpha diversity, and decrease in the relative abundance of the commensal Streptococcus) was observed in esophageal tissues. It should be, however, noted that several bacteriome parameters in patients with BE deviated from this trend. We have also identified some bacterial genera that might (pending further studies) be potentially useful as non-invasive biomarkers or predictors of progression of GERD complications. In particular, periodontal pathogens (such as Porphyromonas) in the oral swabs show promise in this respect, warranting the need for further research in this direction. Very low relative abundance of the genus Faecalibacterium in rectal swabs of patients with EAC may, pending further studies, also potentially serve as an additional indicator of increased risk for EAC development.