The participants were 398 residents (191 men and 207 women) aged 40 years or older, of Shika-machi, Hakui-gun, Ishikawa Prefecture, Japan, whose fecal samples were collected in 2019. The following five conditions were excluded from the analysis. 1) participants without measured serum IgG4 levels, 2) patients taking immunosuppressive drugs as described below; Methotrexate and Enbrel, 3) patients taking medications that significantly affect GM as described below; antibiotics, steroids, bowel regulators and antibacterial agents, proton pump inhibitor (PPI), 4) patients suspected cancer and IgG4-related disease, 5) individuals with missing data, 6) patients with inflammatory bowel disease (IBD).

Data from the Shika-machi Super Preventive Health Examination, a population survey aimed at establishing preventive methods for lifestyle-related diseases, were used. The survey was conducted between 2019. The four model districts selected from the Shika area were Horimatsu, Higashimasuho, Tsuchida, and Higashiki (Karashima et al., 2018; Nagase et al., 2020).

The Shika-machi Super Preventive Health Checkup data regarding parameters such as age, sex, medical history, medication status, allergy status, and alcohol consumption/smoking status were collected using a questionnaire. The body mass index (BMI) was calculated by dividing the current weight (kg) by the square of the height (m²). Venous blood was collected early in the morning after a 12-hour fast. The 24-hour urinary sodium excretion was calculated based on the 24-hour urinary creatinine and sodium excretion values (Nagase et al., 2020). The estimated daily salt intake was calculated using 24-hour urinary sodium excretion.

Immune-related blood samples were measured using the following test kits; IgG4 (IgG4 subclass BS-TIA3 IgG4, MEDICAL & BIOLOGICAL LABORATORIES CO., LTD., Tokyo, Japan); Immunoglobulin G (IgG) (N-assay TIA IgG-SH Nittobo, NITTOBO MEDICAL CO., LTD., Tokyo, Japan); Immunoglobulin E (IgE) (ImmunoCAP Total IgE, THERMO FISHER SCIENTIFIC INC., Waltham, MA, USA); 50% hemolytic unit of complement (CH50) (auto CH50-L eikenII, DENKA COMPANY LIMITED, Tokyo, Japan); Anti-CCP antibody (Stacia MEBLux test CCP, MEDICAL & BIOLOGICAL LABORATORIES CO., LTD., Tokyo, Japan); antinuclear antibody (ANA) (anti-nuclear antibody (ANA) (FA) [FR], FUJIREBIO INC., Tokyo, Japan); Aniti-SS-A/Ro antibody (stacia MEBLux test SS-A, MEDICAL & BIOLOGICAL LABORATORIES CO., LTD., Tokyo, Japan); rheumatoid factor (RF) levels (LZtest ‘eiken’ RF, EIKEN CHEMICAL CO., LTD., Tokyo, Japan).

Fecal samples were collected using previously described methods (Miyajima et al., 2022) and stored at −80°C until DNA extraction. The processing of fecal samples was carried out in a non-proliferation level 2 (P2) laboratory. The DNA extracted from the GM was processed to identify the 16S rRNA gene sequence using a previously reported next-generation sequencing method (Miyajima et al., 2022).

For microbiome analysis, QIIME2 software was used (Bolyen et al., 2022). Demultiplexed paired-end sequence data were denoised with DADA2, and the Silva 16S rRNA database (release 132) naïve Bayes classifier was used for Amplicon Sequence Variant classification (Quast et al., 2013). Samples with fewer than 5000 sequences were removed from the analysis.

Statistical analysis and machine learning were performed using Python (version 3.10.9) (Pedregosa et al., 2011) or R, using R-studio (version 4.2.3, RStudio, Boston, MA, United States).

The clinical information of the participants underwent a normality assessment using the Shapiro-Wilk test. Normally distributed data are expressed as mean ± standard deviation, while non-normally distributed data are presented as median (25^th–75^th percentile). The significance of differences in clinical information between the groups was assessed using Student’s t-test for normally distributed data and the Wilcoxson rank-sum test for non-normally distributed data.

The patients were categorized into two groups, high and low, based on the median values of IgG4. Quade’s non-parametric ANCOVA and Tukey’s HSD test were used to compare the relative proportions of GM between the high and low IgG4 groups. Confounders such as age, sex, BMI, daily salt intake, frequency of alcohol consumption per week, and smoking were adjusted for (Vujkovic-Cvijin et al., 2020). Additionally, the clinical background variables that exhibited significant differences between the high and low IgG4 groups were included as new confounding factors. The significance level for all tests was set at P < 0.05. Alpha diversity was evaluated using the Shannon index, with Amplicon Sequence Variant values (Willis, 2019). To assess the beta diversity, non-metric multidimensional scaling analysis with the Bray-Curtis dissimilarity metric from the “vegan” package in R was used, along with permutation multivariate analysis of variance (Dixon, 2003). To identify GM associated with IgG4, linear discriminant analysis effect size (LEfSe) was employed (Segata et al., 2011).

The odds ratios and P-values were calculated using the least absolute shrinkage and selection operator logistic regression model (LASSO logistic regression) from the “glmnet” package in R (Tibshirani, 1996). Multicollinearity was assessed using the variance expansion factor (VIF) and only bacterial genera with a VIF smaller than 10 were used in the LASSO analysis. Correlation coefficients and P-values were calculated using Spearman’s rank correlation coefficient in R’s “Package ppcor” after adjusting for the variables listed above. The correlation coefficients were plotted using “Package pheatmap”.

The heat maps were visualized as dendrograms using hierarchical clustering, which was based on similarity by correlation coefficient. Bacterial genera that were significantly correlated with IgG4 and one bacterial genus with the closest inter-cluster distance was set up as a new bacterial genus group. The closest bacterial genus was not grouped if it was a population of several bacterial genera.

The direct linear non-Gaussian acyclic model (LiNGAM) model was built using “LiNGAM” in Python (Shimizu et al., 2011). The bacterial genera chosen for LiNGAM algorithm were selected based on their significant associations identified in at least one of the following analyses: ANCOVA, Tukey’s HSD, LEfSe, LASSO, and Correlation analyses. To demonstrate the robustness and consistency of the causal relationships, the occurrence rates and partial regression coefficients of the causal relationships were presented. Unselected bacterial genera were entered exhaustively as noise, and their impact on causality was observed (Mizoguchi et al., 2023).

Data on GM were procured from the fecal specimens of 234 study participants. The study dismissed 138 samples lacking IgG4 quantification, six individuals under immunosuppressive or gut flora-altering medications, one potential cancer case, and one suspected IgG4-RD case. None of the participants had IBD. In total, 88 patients (46 females and 42 males) participated in the analysis. Supplementary Figure 1 contains a flowchart on sample selection. Table 1 elucidates their clinical data. Participants were sorted into high and low categories based on their median serum IgG4 values. The median IgG4 value for all participants was 41.7 mg/dL. By gender, the median values stood at 34.8 mg/dL for females and 57.1 mg/dL for males. Significant discrepancies in BMI, IgG4, IgE, alcohol consumption, and smoking prevalence were observed between genders. Supplementary Tables 1 – 3 respectively offer a comparative overview of the immunological landscape of high and low IgG4 cohorts of all participants, women and men, respectively.

Figure 1A displays stacked plots showing the mean relative abundance of the top 30 bacterial genera among women and men. The top 30 bacterial genera accounted for an average of 81% of women and 85% of men. Figure 1B demonstrates the mean relative abundance of the top 30 genera in the high IgG4 and low IgG4 groups for all participants. Among female participants, the top 30 bacterial genera accounted for an average of 83% of the high IgG4 group and 83% of the low IgG4 group. When segregated by gender, they constituted 83% for both IgG4 groups among women ( Figure 1C ), while for men, they accounted for 82% and 62% in the high and low IgG4 groups, respectively ( Figure 1D ).

Figures 2A–D and 2E–H depict the alpha and beta diversities, respectively, and the analyses revealed no significant disparities in gut GM diversity between sexes ( Figures 2A, E ) or among the high and low IgG4 groups ( Figures 2B–D, F–H ).

Figure 3 demonstrates the significant differences in the presence of specific bacterial genera between IgG4 groups and sexes: Anaerostipes were more prevalent in women than men ( Figure 3A ). In all participants, the proportion of Faecalibacterium present in the High IgG4 group was significantly lower than in the Low IgG4 group, and the proportion of Megasphaera present in the High IgG4 group was significantly higher than in the Low IgG4 group ( Figure 3B ). Women exhibited a significantly lower representation of Anaerostipes in the High IgG4 group ( Figure 3C ), while men in the High IgG4 group had significantly diminished proportions of Faecalibacterium and Ruminococcus.1 ( Figure 3D ).

LEfSe revealed that the relative abundance of Blautia, Megamonas, Prevotella 9, Megasphaera, Ruminococcus.1, Anaerostipes, Subdoligranulum, Escherichia-Shigella has difference among women and men ( Figure 3E ). In all patients, the relative abundance of Megasphaera was higher in the high IgG4 group, and the relative abundances of Anaerostipes and Faecalibacterium were lower in the low IgG4 group ( Figure 3H ). In females, that of Lachnospiraceae was higher in the high IgG4 group ( Figure 3G ). In males, that of Megasphaera was higher in the high IgG4 group, and the relative abundance of Faecalibacterium and [Eubacterium] hallii group was lower in the low IgG4 group ( Figure 3H ).

A predictive model for classifying high/low IgG4 groups was developed using GM data in a LASSO logistic regression model. Thirteen bacterial genera (Blautia, Bifidobacterrium, Subdoligranulum, Streptococcus, Collinsella, Enterobacteriaceae, Fusicatenibacter,. [Eubacterium] hallii group, Anaerostipes, Veillonella, Romboutsia, Lactobacillus, and Ruminococcus.1) had VIF <10 in all participants. In LASSO in all participants, none of the 13 bacterial genera had statistically significant odds ratios. The classification prediction model in all participants was the area under the receiver operating characteristic curve showed 0.658, sensitivity 0.750, and specificity 0.568 ( Figure 4A ).

Twelve bacterial genera (Blautia, Bifidobacterium, Parabacteroides, Collinsella, [Eubacterium] hallii group, Fusicatenibacter, [Eubacterium] coprostanoligenes group, Megasphaera, Anaerostipes, Veillonella, and Romboutsia) had VIF <10 in male. Of the 12 bacterial genera, only Romboutsia was statistically significant with an odds ratio of 2.696 (95% confidence interval 1.031-7.050, P=0.043) ( Supplementary Table 4 ). The classification prediction model in males was the area under the receiver operating characteristic curve showed 0.907, sensitivity 0.857, and specificity 0.857 ( Figure 4B ). Only Parabacteroides had a VIF of less than 10, and other bacterial genera showed multicollinearity in women. For Parabacteroides alone, no predictive model could be built by LASSO and no ROC curve could be drawn.

Figure 5 illustrates the correlation between IgG4 levels and the top 30 intestinal bacterial species, considering the relative abundance. Megasphaera and Lactobacillus displayed significant positive correlations with IgG4 levels in the entire cohort, while Faecalibacterium and [Eubacterium] hallii group exhibited significant negative correlations. In women, IgG4 levels were correlated positively with Megasphaera and negatively with [Eubacterium] hallii group. In men, a negative correlation was observed with Ruminococcus.1.

Based on the similarity of correlation coefficients between IgG4 values and bacterial genera visualized in a dendrogram, three groups were redefined: Group A encompassed Bifidobacterium and Lactobacillus; Group B incorporated Ruminococcus.1 and Ruminococcus.2; and Group C comprised the [Eubacterium] hallii group and Anaerostipes.

Causal inference with the direct LiNGAM model using bacterial genera were significantly associated with IgG4 in ANCOVA, Tukey’s HSD, LefSe, LASSO, and Correlation analyses.

Supplementary Figure 2 shows the results of LiNGAM with bacterial genera and IgG4 listed in Supplementary Table 5 . No causal relationships were estimated between bacterial genera and IgG4.

Figure 6 shows the estimated causal relationship with IgG4, including the redefined bacterial groups. In all participants and women, an increase in Megasphaera was also associated with an increase in serum IgG4 levels, while an increase in group C ([Eubacterium] hallii group and Anaerostipes) was associated with a decrease in serum IgG4 levels. In the robustness analysis of causal results, a causal direction from Megasphaera to IgG4 was detected in 82.6% of cases (partial regression coefficient 385.0 ± 45.9) and from group C to IgG4 in 78.3% of cases (partial regression coefficient -359.4 ± 62.9) in all participants. A causal direction from Megasphaera to IgG4 was detected in 81.8% of cases (partial regression coefficient 673.7 ± 73.7) and from group C to IgG4 in 59.1% of cases (partial regression coefficient -710.0 ± 84.3) in women. In contrast, no causal relationship between bacterial genus and IgG4 could be inferred in males.