The specific design and process of this MR study are divided into 3 steps ( Figure 1 ). In the first and second steps, we utilized bi-directional two-sample univariable Mendelian randomization (UVMR) to assess the causal relationships between exposure and outcome, following the three main assumptions of MR analysis: selected SNPs (1) should be closely associated with exposure (2), affect the outcome only through exposure, and (3) should not be associated with potential confounding factors (10). Reverse analysis reveals results with reverse causal relationships. And validate the results of UVMR using external data. In the third step, we employed both UVMR and multivariable Mendelian randomization (MVMR) to analyze the mediation effects of mediators between exposure and outcome, calculating the effect values and proportions for each qualified mediator. This study adhered to the Strengthening the Reporting of Observational Studies in Epidemiology using Mendelian Randomization (STRBOE-MR) guidelines ( Supplementary Table 1 ) (11).

This study use gut microbiota GWAS data from the Dutch Microbiome Project (DMP) as exposures. This study included 7,738 people of European descent (12). The data was determined through shotgun metagenomic sequencing of fecal samples, encompassing a total of 207 taxonomic units (5 phyla, 10 classes, 13 orders, 26 families, 48 genera, and 105 species).

This study conducted mediation analysis using GWAS summary data of immune cell phenotypes and inflammatory cytokines. The immune cell phenotype data included 3,757 individuals of European ancestry from non-overlapping cohorts, comprising 731 immune features: absolute cell counts (AC, n=118), median fluorescence intensity reflecting surface antigen levels (MFI and SAL, n=389), morphological parameters (MP, n=32), and relative cell counts (RC, n=192) (13). These features, such as MFI, AC, and RC, include mature stages of B cells, CDCs, T cells, monocytes, myeloid cells, TBNK (T cells, B cells, and natural killer cells), and Treg panels. The MP features encompass CDC and TBNK panels. Inflammatory cytokines comprised 41 inflammation-modulating cytokines and CRP. Data for the 41 systemic inflammation-modulating cytokines were obtained from meta-analyses of cytokine-related GWAS from three independent population cohorts, including 8,293 Finnish individuals from the Young Finns Study and the FINRISK studies (FINRISK1997 and FINRISK2002) (14). These 41 cytokines were initially normalized through a first-step inverse transformation of cytokine distributions, followed by a second-step inverse transformation of residuals from linear regression models of transformed cytokines on age, sex, body mass index (BMI), and genetic principal components. CRP data originated from a meta-analysis of GWAS involving 158 European individuals from the Cytokine Working Group (CIWG) consortium, with adjustments made for age, sex, and population structure for 15 genetic associations between 2.4 million genetic variants and log-transformed CRP levels (15).

This study utilized GWAS summary statistics data from the FinnGen Consortium R10 release for four types of pancreatic tumors as outcomes, diagnosis was made using International Classification of Diseases, 10th Revision (ICD-10) diagnosis codes (16). These included Malignant neoplasm of pancreas (MNP, cases=1626, controls=312567) and Benign neoplasm of pancreas (BNP, cases=603, controls=410975), as well as two pathological subtypes: Adenocarcinoma and ductal carcinoma of pancreas (ADCP, cases=731, controls=313462) and Neuroendocrine tumor and carcinoma of pancreas (NTCP, cases=129, controls=314064). Validation of UVMR results External pancreatic cancer data from GWAS Catalog (cases=1196, controls=475049) (17). Please refer to Table 1 for specific details of all the data.

We employed the following steps to select instrumental variables (IVs) for analysis. Firstly, to obtain more effective SNPs, the significance threshold for IVs of GM, IC, and IF was set at 1e-05, while for the IVs of the four pancreatic tumors, it was set at 1e-06. Secondly, PLINK software (version v1.90) was used to remove SNPs with linkage disequilibrium (LD) r2 <0.001 within a distance of 10,000 kb (18). Thirdly, SNPs significantly associated with the outcomes (with a significance threshold of 5e-05) were excluded. Fourthly, palindromic SNPs were removed to ensure that the effect of SNPs on exposure and outcome corresponded to the same allele. Finally, we calculated the F-statistic values to measure the strength of IVs (19), retaining SNPs with F-values greater than 10 and excluding those with a minimum allele frequency (MAF) less than 0.01. Detailed information regarding these IVs can be found in Supplementary Table 2 .

In assessing the causal relationship between GM and the four types of pancreatic tumors, the Inverse Variance Weighted (IVW) method was employed as the primary approach (8), with four other methods used for supplementary analysis (Weighted Median, MR-Egger, Weighted Mode, and Sample mode). When reporting results for binary outcomes, odds ratios (ORs) were presented along with 95% confidence intervals, and when the results were for continuous variables, β values were reported. In the MVMR analysis, the Multivariable Inverse Variance Weighted (MV-IVW) method was used as the primary analysis. Results with a p-value <0.05 were considered significant. Correction for multiple testing of all IVW results was performed using the False Discovery Rate (FDR) method, and FDR q-values were provided. The same approach was used for external validation.

We conducted the screening of potential mediators in the gut microbiota-pancreatic tumor pathway through the following steps ( Figure 2 ). In the first step, the UVMR was used to select the mediators influenced by exposure causally, and their effect values (β1) were calculated. In the second step, the UVMR was used to select the mediators identified in the first step that causally affected the outcome, and their effect values (α) were calculated. In the third step, based on the direction of the obtained exposure-outcome (β), exposure-mediator (β1), and mediator-outcome (α) effect sizes, mediators consistent with logic were retained (if the total effect of exposure on the outcome β is positive, then both β1 and α should be positive or negative; conversely, if the total effect of exposure on the outcome β is negative, then β1 and α should be one positive and one negative). In the fourth step, the MVMR was used to assess the causal effects of mediators on the outcome after adjusting for exposure effects. The mediator with MV-IVW P-value <0.05 is regarded as the final result. Subsequently, combining the causal effect of exposure on the outcome obtained from UVMR (β), the “product of coefficients” method was used to calculate the mediated effect values of mediators in the gut microbiota-pancreatic tumor pathway (β1 × β2) and the proportion of the effect values ([β1 × β2]/β) (20).

We conducted sensitivity analyses using the MR Egger regression method, leave-one-out analysis, and MR-PRESSO method. Cochran’s Q statistic was computed for each SNP to assess heterogeneity, and the p-value from the MR Egger regression intercept test was used to evaluate horizontal pleiotropy (21). The MR-PRESSO method was employed to correct for potential horizontal pleiotropy in the selected IVs (22). P_heterogeneity <0.05 was considered indicative of heterogeneity, while P_Global.test and P_pleiotropy <0.05 were considered indicative of pleiotropy. Evidence of pleiotropy would suggest lack of causal evidence.

All analyses were performed using the TwoSampleMR (23), MR-PRESSO (22), and MendelianRandomization (24) packages in R software (version 4.3.0).

Through the above steps, the number of SNPs for GM ranged from 4 to 60 (median=28), for immune cells ranged from 10 to 189 (median=24), for inflammatory cytokines ranged from 6 to 30 (median=18), and for pancreatic tumors ranged from 201 to 855 (median=622). Additionally, all SNPs had an F-statistic greater than 15, indicating the absence of weak instrumental variables ( Supplementary Table 2 ).

Among the 207 types of gut microbiota included in the analysis, a total of 20 gut microbiota (representing 1 order, 2 families, 2 genera, and 15 species from p_Bacteroidetes, p_Firmicutes, and p_Proteobacteria) were found to have a causal relationship with outcome ( Figure 3A ). Figure 3B presents the results of evaluating the causal effects of gut microbiota on the four types of pancreatic tumors using the IVW method as the primary approach. When ADCP was considered as the outcome, it was observed that s_Ruminococcus_lactaris (OR=0.538; P=0.042) and s_Veillonella_unclassified (OR=0.781; P=0.038) were inversely associated with its risk, while s_Bacteroides_finegoldii (OR=1.318; P=0.009), s_Bacteroides_fragilis (OR=1.325; P=0.018), s_Bacteroides_coprocola (OR=1.464; P=0.008), s_Prevotella_copri (OR=1.529; P=0.031), s_Parabacteroides_distasonis (OR=1.627; P=0.031), o_Pasteurellales (OR=1.649; P=0.035), f_Pasteurellaceae (OR=1.649; P=0.035), and s_Holdemania_unclassified (OR=1.690; P=0.002) were positively associated with its risk. When BNP was considered as the outcome, it was found that s_Bacteroides_finegoldii (OR=1.286; P=0.021), s_Parasutterella_excrementihominis (OR=1.562; P=0.034), and s_Bacteroides_vulgatus (OR=2.003; P=0.011) were positively associated with its risk. When MNP was considered as the outcome, it was found that s_Ruminococcus_lactaris (OR=0.646; P=0.021) was inversely associated with its risk, while s_Bacteroides_ fragilis (OR=1.190; P=0.041), s_Streptococcus _thermophilus (OR=1.285; P=0.040), and s_Holdemania_unclassified (OR=1.418; P=0.003) were positively associated with its risk. When NTCP was considered as the outcome, it was found that s_Paraprevotella_ xylaniphila (OR=0.359; P=0.0004), g_Pseudoflavonifractor (OR=0.372; P=0.038), f_Acidaminococcaceae (OR=0.407; P=0.017), and s_Bacteroides_finegoldii (OR=0.561; P=0.014) were inversely associated with its risk, while s_Holdemania_filiformis (OR=3.369; P=0.019), s_Ruminococcus_torques (OR=4.015; P=0.029), g_Paraprevotella (OR=4.022; P=0.001), and s_Paraprevotella_ unclassified (OR=5.557; P=0.0006) were positively associated with its risk ( Figure 3C ; Supplementary Table 3 ). In addition, we conducted external validation on the gut microbiota that has a causal relationship with pancreatic tumors, Using IVW as the main analysis method, the validation results of the other gut microbiota, except for s_Bacteroidees_vulgatus, s_Prevotella_copri, s_Holdemania_filiformis, and f_Acidamococcaceae, are consistent with the findings. The results are presented in Supplementary Table 3 .

When assessing the reverse causal effects of pancreatic tumors on gut microbiota, using the IVW method as the primary approach, it was found that the abundance of 37 gut microbiota (including 1 phylum, 4 classes, 5 orders, 5 families, 7 genera, and 15 species from p_Actinobacteria, p_Bacteroidetes, p_Firmicutes, p_Proteobacteria, and p_Verrucomicrobia) could be influenced by four types of pancreatic tumors. Specifically, when ADCP was considered as the exposure, it was observed that the abundance of s_Parabacteroides_merdae (β=0.064; P=0.019) could be influenced, showing a contradictory positive causal effect, which violates the fundamental assumption of MR analysis. Therefore, s_Parabacteroides_merdae was not considered to have an effect on ADCP ( Figure 3D ; Supplementary Table 3 ).

Sensitivity analyses were conducted on the results of the above UVMR analyses. All results passed heterogeneity analysis and pleiotropy analysis, with P-values from heterogeneity analysis, MR-Egger intercept, and MR-PRESSO all greater than 0.05, indicating the absence of minimal heterogeneity and horizontal pleiotropy ( Supplementary Table 8 ). Additionally, Leave-one-out analysis showed that in the bidirectional MR analysis, no single SNP significantly altered the causal effects ( Supplementary Figures 1 ,  2 ).

After screening for potential intermediate factors, we identified a total of 1 inflammatory cytokine and 2 immune cell phenotypes that met our selection criteria. Initially, through two-step UVMR and directional screening based on β, β1, and α, we preliminarily selected 9 gut microbiota-immune cell-pancreatic tumor pathways and 1 gut microbiota-inflammatory cytokines-pancreatic tumor pathway, which included 7 immune cell phenotypes and 1 inflammatory cytokine ( Figure 4 ; Supplementary Tables 4 - 7 ). Subsequently, in the fourth step, we used MVMR to evaluate whether the selected intermediates could independently influence the outcomes and calculated the mediation effect values and proportions of the effect values. After adjusting for the influence of gut microbiota, it was found that “SSC-A on HLA DR+ Natural Killer” remained significantly associated with ADCP, “CD4 on naive CD4+ T cell” remained associated with NTCP, and IL-6 remained associated with MNP. Finally, we found that “SSC-A on HLA DR+ Natural Killer” mediated the causal associations between o_Pasteurellales and f_Pasteurellaceae with ADCP, with a mediation proportion of 9.01% (P=0.01), “CD4 on naive CD4+ T cell” mediated the causal association between s_Bacteroides_finegoldii and NTCP, with a mediation proportion of 13.13% (P=0.01), and IL-6 mediated the causal association between s_Streptococcus_thermophilus and MNP, with a mediation proportion of 18.61% (P=0.004) ( Figure 5 ; Table 2 ).

Sensitive analysis was conducted on the UVMR results of the above intermediate analysis. In the analysis of gut microbiota-immune cell pathways, three causal relationships with horizontal pleiotropy were detected and subsequently excluded from our final results. However, the MR-Egger intercept and the P-value of MR-PRESSO for other analysis results were both greater than 0.05, indicating the absence of minimal horizontal pleiotropy ( Supplementary Table 8 ). Heterogeneity analysis revealed potential heterogeneity in some IVs, but the instrument validity test showed that the strength of these IVs was sufficient (all F-statistics ≥ 15; Supplementary Table 2 ), and no horizontal pleiotropy was detected (all P_Global.test and P_pleiotropy > 0.05; Supplementary Table 8 ). Overall, sensitivity analysis confirmed the reliability of these results.