Our MR analysis used gut microbiota features as the exposure data and oral and oropharyngeal cancer data as the outcome ( Figure 1 ). The original investigation was conducted after obtaining ethical approval for each study included in the MR analysis. Genome-wide association study (GWAS) data for the gut microbiota were obtained from the MiBioGen study (Kurilshikov et al., 2021) and the oral and oropharyngeal cancer data were obtained from the MRC IEU OpenGWAS data resource (Elsworth et al., 2020). The original oral and oropharyngeal cancers data were released from OncoArray Oral Cavity and Oropharyngeal Cancer Consortium (Lesseur et al., 2016) and UK Biobank (Burrows et al., 2021). The MiBioGen consortium includes 24 cohorts and totals 18340 participants. In each cohort, the gut microbiome is analyzed through 16S rRNA sequencing using Illumina platforms (MiSeq or HiSeq), while the participants undergo genotyping using whole-genome SNP arrays. All outcomes were based on European ancestry. Each outcome data source we used contained three outcomes: oral cancer, oropharyngeal cancer, and their combination. The OncoArray study comprised a cohort of 6,034 cases and 6,585 controls. The number of cases and controls in North America were 2549 (42.2% of total cases) and 2522 (38.3% of total controls), respectively. In Europe, 2499 cases (41.4% of total cases) and 2928 controls (44.5% of total controls) were reported. For Hispanic or Latin American, OncoArray included 986 cases (16.3%) and 1135 controls (17.2%). Because our MR Study was limited to European populations and gut microbial exposures were not available with respect to Hispanic or Latin American populations, we therefore excluded outcomes if they were Hispanic or Latin American from the study which were restricted to only those of North American and European ancestry. Given the ethnic diversity, Lesseur et al. evaluated associations within continent (Europe, North and South America) using multivariate unconditional logistic regressions under a log-additive genetic model adjusted for age, sex and regional eigenvectors. Genotyping of 13,107 individuals was conducted using Illumina platforms across 12 epidemiological studies. After rigorous quality control measures, a total of 6,034 cases and 6,585 cancer-free controls were included in the final analysis, 4839 cases and 5257 controls of which were of European decent and therefore were used in this MR analysis. Whole-genome imputation was performed using the Haplotype Reference Consortium panel, resulting in approximately 7 million high-quality imputed variants. The UK Biobank began in 2006 by recruiting approximately 500,000 participants between the ages of 38 and 73. Participants completed a series of questionnaires that provided detailed personal and lifestyle information. In addition, participants provided biospecimens, including blood, urine, and saliva, which were subsequently sequenced for genome sequencing and genotyping using the Illumina sequencing platform in the UK Biobank. Details of outcomes are provided in Table 1 . The number of total cases of oral and oropharyngeal cancer combined and controls of the UK Biobank were 839 and 372016, respectively.

We initially extracted instrumental variables using a genome-wide association study p-value threshold of 5e-08. Subsequently, Kurilshikov et al. identified a total of 30 quantitative SNP-feature associations through mbQTL analysis (Kurilshikov et al., 2021). After conducting the main Mendelian randomization (MR) analysis for these exposures, we further performed a secondary analysis using a more lenient threshold criteria, specifically, we extracted five levels of GWAS summary data of exposure (phylum, class, order, family, and genus) derived from MiBioGen, where the number of features were 9, 16, 20, 35, and 131 respectively. Single-nucleotide polymorphisms (SNPs) smaller than the p-value threshold within the locus-wide range (1 × 10^-5) were selected as instrumental variables according to a previous study (Sanna et al., 2019). One MR principle is that there is no linkage disequilibrium (LD) between the included instrumental variables, as a strong LD may lead to biased results, for strong LD may bias the results by genetic confounding or overestimating the association between the genetic variants used and the exposure and outcome. In this study, the SNPs were clumped within each phenotype, clumping processing (R²< 0.1, clumping distance = 500 kb) (Ni et al., 2021) was performed to evaluated LD among the SNPs used to instrument each phenotype separately, where SNPs in LD were removed and the strongest (i.e., that with the smallest p-value) was retained for MR analyses.

The two-sample MR study relies on three critical assumptions to minimize bias. First, the genetic instruments used were robustly associated with exposure. Second, there is no confounding factor between the instrument and the outcome. Third, instrumental variables affected the outcomes only through exposure, meaning that there was no horizontal pleiotropy effect between instrumental variables and outcomes.

When combining summarized data on genetic associations, it is necessary to ensure that genetic associations are expressed per additional copy of the same allele. If a genetic variant is a biallelic single nucleotide polymorphism (SNP) with alleles A and G on the positive strand, then the corresponding base pairs on the negative strand will be T and C. In this case, one dataset may report the association per additional copy of the A allele, and another per additional copy of the T allele – but the same comparison is being made. Allele and strand information can be double-checked by comparing allele frequency information – if the allele frequencies are similar for the A and T alleles, this means that this is a strand mismatch. For palindromic variants – if the alleles were A and T (or C and G), then the same alleles would appear on both the positive and negative strands. In such a case, if the allele frequency is close to 50%, it may be necessary to drop the variant from the analysis if it is not possible to verify that the alleles have been correctly orientated. While this is a conservative recommendation, allele alignment problems have led to incorrect results in MR analyses. Therefore, palindromic SNPs with intermediate (i.e., 0.42%) frequencies were discarded. Analytical methods used included the Wald ratio and inverse variance weighted (IVW) estimators (Burgess et al., 2013). The Wald ratio estimates the causal effect of the exposure on the outcome using a single instrumental variable (IV). When horizontal pleiotropy is absent, IVW results are considered robust. In general, the Wald ratio was used for the analysis with only one SNP and the IVW method was used for the analysis with multiple SNPs The IVW has the strictest assumptions; therefore, this method was considered the main analytical method, to which results from all other methods were compared. Effect estimates and 95% confidence intervals (CIs) reflect the odds ratios (ORs) for disease per standard deviation (SD) higher abundance of microbiome.

In our study, we chose not to apply multiple testing correction in both primary and secondary analyses due to the complex interactions and correlations among microbial features. Although multiple testing correction is common in many studies, we believe that it might be overly stringent in this context, potentially obscuring true causal relationships. This decision is grounded in the complexity of the microbial ecosystem, where microbes often influence each other. Given these interactions, applying multiple testing correction to each microbial feature may be overly conservative and might obscure some signals that have genuine causal relationships.

Colocalization analysis is commonly employed to ascertain whether two phenotypes are driven by the same causal variant within a particular genomic region, thereby bolstering evidence for the association between the two phenotypes. For colocalisation analyses, we extracted SNPs that fell within a 500 kb window upstream and downstream of the SNP used as an instrument for the microbial feature from both the microbiome and oral cancer GWAS datasets. To be more specific, we acquired GWAS summary data from the MRC IEU OpenGWAS data resource. Using the SNPs extracted during the primary MR analysis, we queried their chromosomal number and position within the PheWAS - IEU OpenGWAS project (mrcieu.ac.uk). With this positional information, we adjusted the numerical value by either adding or subtracting 500,000 to establish a positional range spanning 500kb upstream and downstream. With the chromosome number and positional range, we then conducted a filter operation on the GWAS summary data previously downloaded, resulting in the final selection of SNPs located within a 500kb radius both up and downstream of the instrument. The number of extracted SNPs is shown in Supplementary Table 2 .

Specifically, we used the SNP positions that satisfy the threshold criteria and fall within a 500 kb window upstream and downstream as candidate SNPs for extraction in both exposure and outcome traits for colocalization analysis. Colocalization analysis involves five mutually exclusive model assumptions, namely: H0: no association exists between all SNP loci within a genomic region and both the exposure and outcome; H1/H2: significant association exists between the exposure/outcome, respectively, and SNP loci within a genomic region; H3: association exists between both exposure and outcome and SNP loci within a genomic region, driven by distinct causal variants; H4: association exists between both exposure and outcome and SNP loci within a genomic region, driven by the same causal variant. During the colocalization analysis, posterior probabilities (PP.H0-PP.H4) are generated for each of these models. The sum of the posterior probabilities for the five models equals 1. A higher posterior probability for a specific model indicates a higher likelihood of that model assumption being valid given the data. We considered the model assumption to be valid when the posterior probability of that model was greater than 0.80. Priors are as follows: p1: prior probability a SNP is associated with microbiome trait, default is 1e-4; p2: prior probability a SNP is associated with outcome, default is 1e-4; p12: prior probability a SNP is associated with both exposure and outcome, default is 1e-5.

Our colocalization analysis was based on any result obtained from the primary MR analysis that suggested a causal effect of an exposure the OncoArray outcomes. First, we identified SNPs within a 500kb range for both exposure and outcome, and then merged the data based on their rsID, aligning effect alleles. We calculated the variance (VAR) after splitting and cleaning the data by removing missing values, specifying phenotype types: binary phenotype as “cc” (for outcomes) and continuous phenotype as “quant” (for exposure). We conducted colocalization analysis using the “coloc” package and generated locus zoom plots using the “locuscompare” package.

We approximated the variance explained in each gut microbial trait by the respective associated genetic instruments using the following Equation 1 (Shim et al., 2015; Papadimitriou et al., 2020):

In this context, EAF represents the frequency of the effect allele, beta signifies the estimated genetic impact on exposure, N stands for the sample size of the GWAS concerning the association between the SNP and exposure, and SE represents the standard error of the genetic effect. Since EAF was not available for raw data on gut microbes, we inquired the allele frequency (GRCh37) of each SNP through the 1000 genome project (https://www.internationalgenome.org/1000-genomes-browsers/index.html). The power of our MR analyses was assessed using the online calculator mRnd (https://shiny.cnsgenomics.com/mRnd/).

Maximum likelihood (Pierce and Burgess, 2013), MR-Egger (Bowden et al., 2015), weighted median (Bowden et al., 2016), weighted mode, and MR-PRESSO (Verbanck et al., 2018) were used to infer potential causality. The Maximum likelihood (ML) method bears resemblance to IVW, provided there is no presence of heterogeneity and horizontal pleiotropy. If these assumptions hold true, the outcomes will remain unbiased, and the standard errors will be smaller compared to IVW (Yavorska and Burgess, 2017). We used the MR-Egger intercept term to assess pleiotropy. In the presence of horizontal pleiotropy, the MR-Egger effect estimate is a more accurate representation of the unbiased causal effect. If this intercept term was close to zero, then the MR-Egger regression model would be similar to the IVW model. However, if the intercept term was different from zero, it indicates that there may have been horizontal pleiotropy among these IVs. The weighted median can provide consistent estimates of causal effects, providing<50% of SNPs are invalid. We assessed the robustness of the results using sensitivity analysis and tested heterogeneity using Cochran’s Q-test in the IVW test and MR-Egger regression. We conducted MR-PRESSO tests to correct for horizontal pleiotropic effect by removing the IV outliers. MR-Egger regression was used to estimate the effect of pleiotropy; it produced a more robust pleiotropy-corrected causal estimate under the assumptions of no measurement error and instrument strength independent of direct effects (Burgess and Thompson, 2017). Leave-one-out analysis was used to determine whether the original causal effect was driven by a single SNP. We also used the MR Steiger test to examine evidence for directionality in the relationship (i.e., that the exposure caused the outcome rather than the reverse) by comparing the variance explained in the exposure and outcome by the exposure-related IV (Hemani et al., 2017). Furthermore, we assessed instrument strength using the F statistic (Burgess and Thompson, 2011), calculated using the Equation 2:

where R² represents the variance of exposure explained by the selected SNPs, N is the sample size, and k represents the number of instrumental variables. If F< 10, indicating a higher likelihood of weak instrument bias, the association between instrumental variables and exposure was considered weak.

We had three outcome datasets: OncoArray (including North American and European individuals treated as separate data sources) and the UK Biobank. The meta-analysis involved meta-analysing the IVW-specific estimates of the microbial features on the same outcomes across different datasets. The meta-analysis employed both a random-effects model and a fixed-effects model, and was conducted using the meta-package in R. Heterogeneity between studies was assessed using I² and associated p-value. We first determined whether the direction of the effect size for each microbial feature was consistent across the different methods within each data source. If the direction was consistent across methods, we conducted a meta-analysis. If the effect size direction for a particular microbial feature differed among methods for any of the three data sources, then meta-analysis was not performed for that microbial feature.

Analysis was performed using R software (version 4.0.2; R Foundation for Statistical Computing, Vienna, Austria), with the R packages “TwoSampleMR” (version 0.5.6) (Hemani et al., 2018) and “MRPRESSO” (Verbanck et al., 2018). We followed STROBE-MR (guidelines for strengthening the reporting of observational studies in epidemiological studies using MR to report our results (Skrivankova et al., 2021).

We employed a genome-wide association study p-value threshold of 5e-08 to extract instrumental variables. Our results are consistent with those of Kurilshikov et al., as a total of 30 SNPs associated with 27 microbial phenotypes were extracted (Kurilshikov et al., 2021). Among these, genus Intestinibacter, genus RuminococcaceaeUCG009, and genus CandidatusSoleaferrea yielded one extracted SNP each (rs10805326, rs8009993, and rs830151, respectively), yet no corresponding SNP was found in the outcome. Conversely, for genus Peptococcus, genus Enteror-habdus, genus Ruminococcus1, and genus Faecalibacterium, the sole SNPs (rs75754569, rs11098863, rs10769159, and rs12320842, respectively) were discarded for being palindromic with intermediate allele frequencies.

We proceeded to perform MR analysis on the remaining 20 gut microbial features. Wald ratio results indicated that order Gastranaerophilales, family Gastranaerophilales, and genus Gastranaerophilales exhibited carcinogenic effects on oral and oropharyngeal cancer in the North American population from OncoArray with odds ratios (ORs) per standard deviation (SD) higher abundance of 2.610 [1.124-6.061] (oral and oropharyngeal cancer) and 3.609 [1.232-10.568] (oropharyngeal cancer). However, these effects were inconsistent across multiple data sources, as the causal effect directions of order Gastranaerophilales, family Gastranaerophilales, and genus Gastranaerophilales varied across different data sources ( Supplementary Table 1 ). Subsequently, we conducted colocalization analysis on these 30 SNPs, revealing that none of the PP.H4 values exceeded 0.80, providing limited evidence in support of the hypothesis that “exposure and outcome are associated with SNP loci within a genomic region, driven by the same causal variant.”( Supplementary Table 2 ) ( Supplementary Figure 1 ).

Subsequently, we conducted analyses using a more lenient p-value threshold. We extracted five levels of GWAS summary data on exposure derived from MiBioGen. After the removal of palindromic SNPs, clumping, and harmonization of the data, the number of SNPs associated with the microbial traits that were available in the outcome data ranged from 4 to 22; the microbial trait associated with the fewest SNPs was the genus Hungatella, with four SNPs extracted, and the microbial trait associated with the most SNPs was the class Actinobacteria, with 22 SNPs extracted ( Supplementary Table 3 ). The R² and F-statistic values for the exposures are summarized in Supplementary Table 4 . The R² values for all SNPs of exposure were greater than the corresponding R² for the SNPs of the outcome. The F-statistics for all exposures were greater than 10, indicating little evidence for weak instrument bias, which demonstrated that all SNPs had sufficient validity.

A secondary MR analysis provided evidence for causal relationships between 90 (includes 58 microbiota signatures) and oral cancers across the three data sources. There were 13, 15, and 17 microbiota features that caused a change in risk of oral cancer, oropharyngeal cancer, and combined oral and oropharyngeal cancer in the European population, respectively ( Supplementary Figures 2-13 ), and nine, three, and eight microbiota features that caused a change in risk of oral cancer, oropharyngeal cancer, and combined oral and oropharyngeal cancer in the North American population, respectively ( Supplementary Figures 14-25 ). The number of features that possessed causal relationships with oral cancer, oropharyngeal cancer, and combined oral and oropharyngeal cancer in the UK Biobank data were 10, 5, and 10, respectively ( Supplementary Figures 26-37 ) ( Supplementary Table 3 ).

We present the results of microbiota signatures with causal effects on oral or oropharyngeal cancer in combined oral and oropharyngeal cancer data, for these estimates were consistent in terms of direction of effect in all methods. For the European region, genus Hungatella (OR=2.735, [1.565;4.78] (Oropharyngeal cancer)), and Parabacteroides were risk factors, whereas the genus Alistipes and Eubacterium coprostanoligenes group, phylum Lentisphaerae (OR=0.448, [0.292;0.686] (Oropharyngeal cancer)), class Lentisphaeria, and order Victivallales were protective factors ( Figure 2A ).

In North American populations, one protective factor (genus Butyrivibrio) (OR=0.784, [0.619;0.991] (Oral cavity and oropharyngeal cancer)), and one risk factor (genus Eubacterium xylanophilum group) (OR=2.223, [1.129;4.376] (Oral cavity cancer)) were causally associated with combined oral and oropharyngeal cancers and oral or oropharyngeal cancer ( Figure 2A ).

In the UK Biobank data, there were seven microbial characteristics were associated with outcomes, among which class Betaproteobacteria, order Burkholderiales(OR=1.001, [1.000;1.002] (Oral cavity cancer)) were risk factors, whilst class Actinobacteria, family Bifidobacteriaceae, genus Bifidobacterium(OR=0.998, [0.998;0.999](Oral and oropharyngeal cancer)) and Butyrivibrio, and order Bifidobacteriales were protective factors. ( Figure 2B ).

Unfortunately, none of the signatures were consistent across the three data sources. However, there were 7 signatures that were consistent across the two data sources, including class Lentisphaeria (European region and North American region of OncoArray), genus Butyrivibrio (North American region of OncoArray and UK Biobank), order Bacillales (European region of OncoArray and UK Biobank), order Bifidobacteriales (European region of OncoArray and UK Biobank), order Burkholderiales (North American region of OncoArray and UK Biobank), order Victivallales (European region and North American region of OncoArray), family Bifidobacteriaceae (European region of OncoArray and UK Biobank). We compared the direction of the resulting effect size for the 20 exposures in the main analysis with the direction of the effect size obtained after the use of the lenient threshold criteria. Among the 180 causal relationships examined, 113 exhibited consistent effect directions in the secondary MR analysis following lenient threshold criteria. Of the 20 phenotypes assessed in the main analysis, none had directionally consistent effect estimates across main and secondary analyses across all methods in 9 outcomes. Notably, genus Tyzzerella3 and genus Romboutsia maintained consistent directions for 8 effect estimates after the secondary analysis. The direction of effect for Gastranaerophilales on oral cavity and pharyngeal cancer and oropharyngeal cancer (Both from North American region of OncoArray) was the same as the direction of effect for Gastranaerophilales in the secondary analysis across all methods ( Supplementary Table 1 ).

The results of Cochran’s Q test suggested that there was little evidence of heterogeneity, where all p-values were greater than 0.05. The horizontal pleiotropy between instrumental variables and outcomes was assessed using the MR-Egger regression. There was little evidence of horizontal pleiotropy was found ( Supplementary Table 4 ). MR-PRESSO global test results showed the presence of outliers only in the genus Veillonella analysis of Oral and oropharyngeal cancer (UK Biobank). After excluding outliers, the MR-PRESSO results showed that there was no evidence for a causal effect of Veillonella on either cancer. Additionally, the leave-one-out analysis showed that none of the identified causal associations were driven by any single IV. ( Supplementary Figures 2–37 ).

In our secondary MR Analysis, 18 causal relationships (17 microbial traits) were found to be inconsistent in the direction of causal effects among different methodologies ( Supplementary Table 3 ), so these causal relationships were excluded from the subsequent meta-analysis. There were 72 causal relationships (50 microbial traits) in the same direction of causal effects across methodologies ( Supplementary Table 3 ).

We performed the MR Steiger directionality test to determine whether there was evidence that the assumed direction of causality was correct. All Steiger test p-values were less than 0.05, suggesting that the assumed direction of effect from the gut microbiome to oral and oropharyngeal cancers was correct, except for the genus Hungatella on oropharyngeal cancer (Steiger P-value=0.32). However, the R² of exposure SNPs was greater than that of the outcome, suggesting that all causal directions were correct.

For class Lentisphaeria, genus Butyrivibrio, genus Ruminococcus 2, order Victivallales, and phylum Lentisphaerae consistently aligned with the MR methods, and therefore, estimates of the effect of these microbial features on oral and oropharyngeal cancer from all three data sources were included in the meta-analysis. Our results showed that there was limited evidence suggesting a causal effect of any microbial feature on oral and oropharyngeal cancer in the meta-analysis ( Supplementary Table 5 ). Next, we conducted a meta-analysis for oral cancer outcomes, using the same criteria as before. Although genus Eubacterium xylanophilum group, genus Parabacteroides, and genus Ruminococcus gauvreauii group consistently aligned with the MR methods, the meta-analysis results provided little evidence for causal effects ( Supplementary Table 6 ).

Only the estimates of genus Ruminococcaceae UCG010 were consistent in direction across all methods in individual data sources, therefore, the meta-analysis was conducted with this one microbial feature. Results showed little evidence for a causal effect of Ruminococcaceae UCG010 on oropharyngeal cancer ( Supplementary Table 7 ).