This is a two-sample MR study that rests on three key assumptions for inferring causality (Smith and Ebrahim, 2003; Lawlor et al., 2008). First, the genetic instruments strongly predict the exposure of interest (i.e., the relevance assumption). Second, the genetic instruments are independent of confounders of the exposure–outcome association (i.e., the independence assumption). Third, the genetic instruments affect the outcome only via the exposure of interest (i.e., the exclusion–restriction assumption).

We extracted genetic instruments, i.e., single-nucleotide polymorphism (SNP), for 21 genera of gut microbiota from summary statistics in MiBioGen, which is the largest and latest available 16S fecal microbiome data (up to 18,340 participants, including 16,632 adolescents and 1,708 children) on gut microbiota of individuals dominated by of European descent (∼72.3%, 13,266 participants) (Kurilshikov et al., 2021). We selected instruments suggestively [i.e., p < 5 × 10^–6, which is commonly used to highlight “suggestive” genetic variants (Manolio, 2010)] and independently (i.e., r² < 0.01, based on the 1000G European population reference panel) associated with each genus of gut microbiota. The median age of participants in MiBioGen was ∼46.5 years (range: 4–89 years) with ∼44.2% (∼8,073 participants) of males. The genetic effects on genera of gut microbiota were meta-analyzed using either the abundance levels (mbQTL, i.e., samples with zero abundance were truncated) or a binary trait coding presence/absence (mbBTL, i.e., presence versus absence of the bacterial genus), depending on their absolute abundance levels across sub-cohorts. These microbiome GWASs adjusted for age, sex, study-specific covariates, and the top genetic principal components for population stratification. Additional file 1: Supplementary Table S1 shows the list of genetic instruments for genera of gut microbiota included in this study.

We extracted summary statistics of BE from the most recent meta-analyzed BE GWAS (access via study accession GCST90000515 (Ong et al., 2021), including 3,513 cases and 14,052 controls in UK Biobank (An et al., 2019); and 6,167 cases and 17,159 controls in a sub-meta BE GWAS (Gharahkhani et al., 2016)) based on the rs number of the identified genetic instruments for 21 genera. Specifically, for UK Biobank, BE and EAC diagnosed with ICD-10 codes (i.e., K22.7 for BE and C15 for EAC) based on self-report and clinical diagnosis were combined as one phenotype mainly because BE is the premalignant precursor of EAC and has a very high genetic correlation with EAC (An et al., 2019). For the sub-meta BE GWAS of all individuals of European descent, patients with BE were identified by the histopathological diagnosis of intestinal metaplasia (Gharahkhani et al., 2016). The meta-analyzed BE GWAS adjusted for age, sex, study-specific covariates, and the first 10 principal components. Additional file 1: Supplementary Table S1 shows the genetic instrument associations with BE.

Given that previous studies (including MR and systematic reviews) showed potential roles of obesity (Thrift et al., 2014), depression (Ong et al., 2021), years of schooling (Ong et al., 2021), and cigarette smoking (Andrici et al., 2013; Cook et al., 2012), we further explored associations of the genetic instruments strongly predicted genera of gut microbiota in respective GWASs with obesity proxied by the waist–hip ratio (GIANT and UK Biobank, n = 697,924 participants) (Pulit et al., 2019), cigarette smoking proxied by cigarettes smoked per day (a meta-analyzed GWAS, n = 337,334 participants) (Liu et al., 2019), depression (a meta-analyzed GWAS, n = 246,363 cases and 561,191 controls) (Howard et al., 2019), and years of schooling (SSGAC, up to 1,131,881 participants) (Lee et al., 2018). Details about study participants included in these GWASs are presented in Additional file 1: Supplementary Table S2. We excluded the instrument associated with any of these phenotypes aforementioned at genome-wide significance (p < 5 × 10⁻⁸) to reduce the possibility of pleiotropy. Here, we did not consider the potential pleiotropic effects of GERD on BE because they shared considerable genetic variants, responding for ∼91% (An et al., 2019). In addition, we cross-checked the associations of the selected instruments with any causes of death using the comprehensively genotype-to-phenotype cross-reference PhenoScanner at p < 5 × 10⁻⁸ to reduce other sources of pleiotropy and the risk of selection bias (Kamat et al., 2019; Schooling et al., 2020; Yang et al., 2021). Additional file 1: Supplementary Table S3 shows the excluded instruments and the corresponding pleiotropic effects.

We assessed the instrument strength using the F statistic, whose value less than 10 indicated a higher likelihood of weak instrument bias (Bowden et al., 2016a). We used the multiplicative random-effects inverse-variance weighted method to estimate the association of genetically predicted higher genus of gut microbiota on BE by pooling the Wald estimator [i.e., the ratio between the SNP-outcome effect and the SNP-exposure effect with its standard error being approximated using the first-order weights (Bowden et al., 2017)] of each SNP. We reported the heterogeneity of the Wald estimator using the Cochrane Q statistics (Bowden et al., 2017), with the potential directional pleiotropy indicated by p<0.05 for MR-Egger intercept (Bowden et al., 2015).

We also conducted the sensitivity analysis using a weighted median estimator (Bowden et al., 2016b), MR-Egger (Bowden et al., 2015; Burgess and Thompson, 2017), MR-PRESSO (Verbanck et al., 2018), and MR robust adjusted profile score (MR-RAPS) (Zhao et al., 2020). Specifically, the weighted median method allows for up to 50% of weights from invalid instruments and produces a consistent causal estimate (Bowden et al., 2016b). MR-Egger detects the directional pleiotropy using the p value for the intercept and provides causal estimates after adjusting for pleiotropic effects with an additional assumption about the instrument strength independent of the dependent effect (i.e., InSIDE assumption) (Bowden et al., 2015; Burgess and Thompson, 2017). The MR-PRESSO detects the outlying instruments and provides consistent causal estimates after removing possible outliers with an additional assumption of no invalid instruments (Verbanck et al., 2018). The MR-RAPS allows the inclusion of weak instruments and the presence of systematic and idiosyncratic pleiotropy and provides consistent causal estimates using an adjusted profile likelihood estimator (Zhao et al., 2020).

We approximated the variance of each genus of gut microbiota explained by the included instruments using well-established methods for continuous and binary exposures (Lee et al., 2012; Yarmolinsky et al., 2018). We assessed the power of our MR analyses using the online calculator mRnd (https://shiny.cnsgenomics.com/mRnd/) (Brion et al., 2013). Additional file 1: Supplementary Table S1 presents the estimated variance of each genus explained by the instruments. Additional file: Supplementary Table S4 shows the power of the estimated gut microbiota–BE associations.

All analyses were performed using R Version 3.6.2 (R Core Team (2019). R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/), with the R package “TwoSampleMR” (Hemani et al., 2018). We reported a two-sided p value at the Bonferroni-corrected threshold of 0.05/22 = 0.002, and the p value between 0.002 and 0.05 was considered suggestive of causation. We adhere to STROBE-MR: Guidelines for strengthening the reporting of observational studies in epidemiological studies using Mendelian randomization for reporting our results (Skrivankova et al., 2021).

This analysis of publicly available data does not require ethical approval.

This MR study, taking advantage of the largest publicly available GWAS on gut microbiota and BE, was the first extensive analysis investigating the potential roles of a broad range of the dominant genera of human gut microbiota in BE. Our study found that genetically predicted higher Actinomyces and higher Ruminiclostridium appeared to protect against BE. Our study also found that genetically predicted higher Alistipes, higher Eubacterium, and higher Veillonella were suggestively associated with lower BE risk, while higher Faecalibacterium was associated with higher BE risk.

In particular, there was no evidence supporting an association of Streptococcus with BE in our study, in which Streptococcus was thought to be the dominant BE-specific microbiota (Yang et al., 2009; Gall et al., 2015). This finding is inconsistent with that of previous observational studies showing that higher Streptococcus was positively associated with BE (Lopetuso et al., 2020). However, the discrepancies between previous findings and our study may be due to the chance finding caused by the small sample size and the unmeasured confounding caused by other factors, such as obesity and health status.

Our findings were consistent with those of previous studies showing that higher Actinomyces and higher Veillonella were inversely associated with BE (Elliott et al., 2017; Snider et al., 2019; Zhou et al., 2020), although the Veillonella-BE association became less evident after adjusting for multiplicity. However, our findings did not show any association of Prevotella with BE risk, which was inconsistent with previous observational studies showing a possible association of Prevotella with BE/EAC (Elliott et al., 2017; Snider et al., 2019; Lopetuso et al., 2020). Such discrepancies may be the consequence of confounding bias in observational studies arising from the diet, environmental factors, or medication (Rothschild et al., 2018; Hills et al., 2019), which were thought to shape the composition and abundance of genera of gut microbiota in humans in real-time, probably through mediations of metabolites and inflammatory cytokines (e.g., IL-8) (Munch et al., 2019).

Our findings also suggested that higher Alistipes appeared to protect against BE, which was not reported in previous studies. Nevertheless, from a biological perspective, higher Alistipes (i.e., a recently discovered gram-negative and anaerobic genus of the Bacteroidetes phylum in mostly healthy human GI tract) acts as a potential pathogen, contributing to the beneficial immunomodulation in cancer [e.g., colorectal cancer (Moschen et al., 2016)]. Second, Alistipes has been reported to modulate the tumor microenvironment and gut inflammation (Parker et al., 2020), and thus might have a role in cancer immunotherapy. For instance, manipulating the tumor microenvironment by reducing tumor necrosis factor produced myeloid cells using antibiotics appeared to decrease the tumor eradication rate (Iida et al., 2013). Furthermore, patients with non–small-cell lung cancer who responded to nivolumab (a checkpoint inhibitor for PD-1) tended to have an elevated Alistipes (Jin et al., 2019). Third, Alistipes-related dysbiosis has already been implicated in several other diseases (e.g., liver fibrosis (Rau et al., 2018), cardiovascular diseases (Zuo et al., 2019), and mood disorder (Bangsgaard Bendtsen et al., 2012)), possibly via inflammation.

In addition, our findings suggested that higher Faecalibacterium [i.e., a gram-positive and anaerobic bacterium that is one of the most abundant and critical commensal bacteria of human gut microbiota in the intestine (Miquel et al., 2013; Bag et al., 2017)] appeared to increase BE, possibly due to its role in boosting the immune system (Miquel et al., 2013) and improving gut barrier function (Stenman et al., 2016). Similar effects were also reported in Crohn’s disease (Wright et al., 2015). Finally, our findings suggested that higher Eubacterium (a gram-positive genus of the Eubacteriaceae family) and higher Ruminiclostridium may also have roles in reducing BE risk. However, further research studies using randomized controlled trials remain required to verify such findings.

Although MR provides less confounded estimates of the gut microbiota–BE associations, limitations still exist. First, MR estimates rest on stringent assumptions (e.g., the independence and exclusion–restriction assumptions), which are always untestable (Glymour et al., 2012). However, in this study, we selected genetic instruments that strongly predicted gut microbiota from the largest publicly available GWAS, with an F statistic being greater than 10 to reduce the possibility of weak instrument bias (Bowden et al., 2016a). Furthermore, we selected valid genetic instruments at a “suggestive” threshold of p < 5 × 10^–6, instead of the traditional genome-wide significance (i.e., p < 5 × 10^–8), which may induce weak instrument bias. However, the sensitivity analyses, particularly MR-RAPS allowing the inclusion of weak genetic instruments (Zhao et al., 2020), yielded consistent results, indicating a less likely weak instrument bias. We excluded instruments associated with any potential confounders of the gut microbiota–BE associations and any known pleiotropic effects to reduce the risk of violating the independence and exclusion–restriction assumptions. We also conducted Cochran’s Q statistic test to detect the potential heterogeneity of the causal estimates and MR-Egger regression to examine the possible pleiotropic effects. No heterogeneity and pleiotropy were noted, as shown in Figures 1–3 and Additional file 1: Supplementary Table S4.

Second, the effects of genetic instruments on gut microbiota composition and abundance vary considerably across GWAS studies, and seldom do replications across these studies be available (Kurilshikov et al., 2017). As such, population-specific microbiota compositions may exist, which may result in a broad uncertainty and a lack of reproduction of our findings. However, the use of the largest publicly available GWAS of MiBioGen and extensive sensitivity analyses may reduce the risk of such issues. Furthermore, genetic instruments typically explain a small variation in a specific genus of gut microbiota, inducing underpowered MR estimates. However, our power analysis shows adequate power (i.e., 80%) to detect a positive association at the α = 0.05 level, given the true association exists.

Third, we limited MR analyses to participants of mainly European ancestry to reduce the risk of population stratification. Thus, our findings may not extend to other populations, although causations are always consistent across populations (Lopez et al., 2019). In addition, summary statistics obtained from GWAS on genera of gut microbiota had partially overlapping sets of participants from UK Biobank, which may bias our estimates with uncertain magnitudes and directions (Burgess et al., 2016). However, such an impact seemed small as MR-RAPS robust to the overlapping sample issue yielded similar results (Zhao et al., 2020).

Fourth, canalization buffering genetic factors may also exist. However, its impact on our estimates remains unknown. Fifth, our causal estimates reflected the natural genetic variation in lifelong exposure of genera of gut microbiota on BE, which might be different from the effects of short-term interventions with antibiotics.

Nevertheless, our findings provide genetic evidence for the associations of several genera of human gut microbiota with the risk of developing BE. Replication using different populations when the large-scale GWAS becomes available remains valuable to cross-validate our findings. Furthermore, manipulation of the abundance of Actinomyces and Ruminiclostridium may be helpful in preventing BE. However, Actinomyces are opportunistic pathogens in humans and may cause the infection of actinomycosis (Kononen and Wade, 2015), which should be considered. In addition, future research involving the associations of gut microbiota and BE might prioritize these four genera of gut microbiota (i.e., Alistipes, Eubacterium, Veillonella, and Faecalibacterium).