The overall design of this study is shown in Figure 1 . Mendelian randomisation studies must meet three core assumptions of correlation, independence and exclusivity, namely: (1) instrumental variables must be highly correlated with exposure factors, (2) instrumental variables must not be correlated with any confounding factors associated with ‘exposure-outcome’, and (3) instrumental variables can only influence outcome variables through exposure factors (Burgess et al., 2017).

Genetic association data for the gut microbiota were obtained from a large GWAS study by the MiBioGen consortium (https://mibiogen.gcc.rug.nl), which included 18,340 participants from 11 countries. After removing 15 unknown bacterial taxa, we ended up with 9 phyla, 16 classes, 20 orders, 32 families and 119 genera, for a total of 196 taxa ( Table S1 ) (Kurilshikov et al., 2021). GWAS summary statistics for BD were obtained from the IEU GWAS database (https://gwas.mrcieu.ac.uk/), a dataset (ieu-b-41) containing 20,352 cases of European ancestry and 31,358 controls of European ancestry. All patients with BD fulfilled international consensus criteria (DSM-IV or ICD-10) (Stahl et al., 2019).

First, we extracted SNPs closely associated with gut microbiota from the GWAS study data using a threshold of P<1×10^-5 (Zeng et al., 2023). In addition, to ensure the independence of each IV, we used a threshold of r²<0.001 and a window size of 10,000 kb to mitigate linkage disequilibrium (Zeng et al., 2023). At the same time, palindromic SNPs and SNPs not present in the results were also removed from the IVs. Finally, SNPs with an F-statistic <10 were removed to eliminate bias caused by weak instrumental variables in the results (Burgess and Thompson, 2011). F=β² (exposure)/SE² (exposure) was used to calculate the strength of the IV (Lawlor et al., 2008; Burgess et al., 2017; Zhang et al., 2022).

Inverse-variance-weighted (IVW) method, weighted median method, weighted model, MR-Egger and Mendelian randomization pleiotropy residual sum and outlier (MR-PRESSO) test were used in this study for Mendelian randomization analysis. The IVW method uses Wald estimators and delta’s to compute ratio estimates for each SNP, and then combines estimates from each SNP to obtain primary causation estimates (Burgess et al., 2013). Compared with the other four methods, IVW provided more accurate effect estimates and was therefore chosen as the primary method of analysis, with the other methods used as supplementary validation of IVW (Bowden et al., 2016; Larsson et al., 2017; Yavorska and Burgess, 2017). The MR online power calculation tool (https://shiny.cnsgenomics.com/mRnd/) was used to calculate the statistical power of the causal effect estimates (Burgess, 2014).

In addition, as the Inverse Variance Weighted method may be subject to invalid instrument bias or pleiotropy, this study tested the validity and robustness of the IVW results through sensitivity analysis. Firstly, the potential heterogeneity of the MR analysis results was quantified and tested using Cochran’s Q. Secondly, the horizontal pleiotropy of the results was assessed using the MR Egger intercept test and the MR PRESSO global test, with p < 0.05 as statistically significant (Burgess and Thompson, 2017). By setting the number of distributions to 10,000, the MR-PRESSO outlier test was also used to adjust for horizontal pleiotropy by detecting and removing outliers (Verbanck et al., 2018). In parallel, we assessed the effect of abnormal SNPs on the results of Mendelian randomization using leave-one-out analysis. To further rule out confounding, instrumental variables with significant MR estimates related to risk factors for bipolar disorder (smoking, obesity, type 2 diabetes) were examined and excluded using the PhenoScanner, and re-tested to see if the causal effects remained significant (Kamat et al., 2019; Vermeulen et al., 2021; Miola et al., 2022).

To obtain more rigorous conclusions, we corrected the P-values using the Bonferroni method. Considering that each character level (phylum, class, order, family and genus) includes several gut microbial taxa, the corrected threshold for each character level should be 0.05/N, where N is the number of taxa included in that character level. Thus, the corrected thresholds for the phylum, class, order, family and genus levels are 5.56 × 10^-3 (0.05/9), 3.13 × 10^-3 (0.05/16), 2.50 × 10^-3 (0.05/20), 1.56 × 10^-3 (0.05/32) and 4.20 × 10^-4 (0.05/119) respectively. MR results with p-values less than the Bonferroni correction threshold were considered significant. MR results with p-values < 0.05 were considered nominally significant (Liu et al., 2023; Luo et al., 2023). A p-value < 0.05 was considered statistically significant for sensitivity analysis in this study. The association between gut microbiota and bipolar disorder was expressed as an odds ratio (OR) with its 95% confidence interval (CI).

All of the above analyses were primarily performed using the Two-Sample-MR package (version 0.5.7) with R software (version 4.2.3).

Through rigorous screening in the previous section, 2527 SNPs were identified as instrumental variables. The phylum, class, order, family and genus levels of the gut microbiota contained 122, 221, 272, 432 and 1480 instrumental variables, respectively. Details of the IVs are provided in Table S2 .

The results of the preliminary analysis of the relationship between gut microbiota and bipolar disorder are shown in Figure 2 and Table S3 .

We used IVW as the primary method for MR analysis of the 196 taxa tested. We detected a total of 16 taxa with P values < 0.05 for MR analysis results, and of these 16 taxa with nominal significance, we did not find any taxa with P values less than the Bonferroni correction threshold for MR analysis results. As shown in the sensitivity analysis below, the Ruminococcaceae NK4A214 group was found to be horizontally pleiotropic in the MR-Egger intercept test. The Ruminococcaceae UCG-003 and Bacillales were found to be horizontally pleiotropic in the MR-PRESSO global test. These three taxa are therefore removed from the 16 taxa of nominal significance. So after removing these 3 taxa, we ended up with 13 nominally significant gut microbiota. Prevotella 7 (OR=1.09, 95%CI=1.01-1.17, P=0.027), Peptococcus (OR=1.12, 95%CI=1.04-1.21, P=0.003), Butyricimonas (OR=1.15, 95%CI=1.02-1.29, P=0.021), Desulfovibrionales (OR=1.15, 95%CI=1.00-1.32, P=0.049), Eubacterium xylanophilum group (OR=1.16, 95%CI=1.01-1.33, P=0.032), Betaproteobacteria (OR=1.17, 95%CI=1.02-1.33, P=0.027), Acidaminococcaceae (OR=1.19, 95%CI=1.02-1.40, P=0.026), Roseburia (OR=1.20, 95%CI=1.03-1.41, P=0.018), Burkholderiales (OR=1.22, 95%CI=1.05-1.41, P=0.008) and Terrisporobacter (OR=1.25, 95%CI=1.04-1.49, P=0.016) had a potential causal effect on the increased risk of BD, whereas Ruminiclostridium 5 (OR=0.86, 95%CI=0.75-0.98, P=0.028), Candidatus Soleaferrea (OR=0.91, 95%CI=0.84-0.98, P=0.020), Victivallis (OR=0.92, 95% CI=0.86-1.00, P=0.036) had a potential causal effect on the decreased risk of BD ( Figure 3 ). The MR results of the Weighted Median, Weighted Mode, MR-Egger and MR-PRESSO methods were then used as additional validation of the IVW. Of the 13 taxa with solid MR results in IVW, four taxa (Betaproteobacteria, Butyricimonas, Candidatus Soleaferrea and Burkholderiales) showed MR-Egger results that were not parallel to those in IVW. Normally, in this case, we would need to re-examine the instrumental variables at more stringent thresholds (Chen et al., 2021). However, these four taxa already contain few instrumental variables, and it would be difficult to ensure the statistical power of the MR results if the number of instrumental variables were further reduced. 1 × 10^-5 is the threshold chosen for most Mendelian randomization studies of the gut microbiota (Yu et al., 2023; Zeng et al., 2023). In addition, the causal effects of the IVW analysis were more precise than those of the MR-Egger analysis. Finally, this situation was already observed in the study by Luo et al. who ultimately chose to retain taxa for which the MR-Egger results did not parallel the IVW results (Luo et al., 2023). Thus, in the absence of heterogeneity and horizontal pleiotropy, it was ultimately acceptable to retain the IVW results for these four taxa.

As MR results can be affected by invalid instrument bias or pleiotropy, this study examines the validity and robustness of the results through sensitivity analysis. With the exception of the Ruminococcaceae NK4A214 group, Ruminococcaceae UCG-003 and Bacillales, none of the 13 gut microbial taxa were found to be horizontally pleiotropic in the MR-Egger intercept test and the MR-PRESSO global test ( Table 1 ). Figure 4 visually demonstrates that the MR results are not affected by horizontal pleiotropy. For these 13 taxa, the Cochran’s Q test and funnel plots provided no evidence for the presence of heterogeneity ( Figure S1 ). In addition, the leave-one-out analysis demonstrates the robustness of the MR results, as the exclusion of any of the instrumental variables does not fundamentally affect the results ( Figure S2 ). The visualisation results of the remaining analyses are shown in Figures S3, 4 .

To further rule out confounding, we examined instrumental variables of 13 gut microbiota taxa associated with bipolar disorder using the PhenoScanner. Of these, rs6058181, rs113054641, rs72814525, rs12500231, rs56349194 were associated with obesity, rs45497800 with smoking and rs6494306 with type 2 diabetes. After removing these instrumental variables associated with confounders, the causal relationships between the 13 gut microbial taxa and BD were reassessed. The results showed that the causal relationships of the remaining gut microbial taxa remained significant, with the exception of Acidaminococcaceae and Victivallis ( Table S4 ; Figure 5 ).