Participants in this study (Takeuchi et al., 2023) were recruited from 2014 to 2016 during their annual health check-ups at Tokyo University Hospital. The participants were Japanese individuals aged between 20 and 75 years, including both males and females. Exclusion criteria was applied, such as a prior diagnosis of diabetes, routine use of diabetes or intestinal medications, recent antibiotic use, and significant weight loss in the three months prior to sample collection. To ensure comparable clinical characteristics, the study enrolled 101 individuals with normal health, 100 individuals classified as obese (based on the Japanese definition of Body Mass Index [BMI] ≥ 25), and 112 individuals classified as prediabetic (based on FBG ≥ 110 mg/dL and/or HbA1c ≥ 6.0%) using their clinical data. Participants were instructed to fast overnight before their hospital visit, during which clinical information and blood samples were taken in the morning. Blood samples were immediately processed and stored at −80°C. Fecal samples were also collected in the morning, transported to the hospital within 24 hours, and stored at −80°C. Out of the total, 256 participants provided fecal samples on the day of their hospital visit, while the rest collected their samples between 2 days before and 7 days after the visit. A small number of participants either collected their samples too long after the visit, collected the evening before the visit, or did not provide fecal samples. Two individuals withdrew from the study after enrollment. Therefore, a total of 306 individuals underwent physical examination, laboratory tests, and fecal sampling. Fecal metagenomic data were available for 290 individuals due to limited samples.

Genomic DNA was extracted from peripheral blood samples of 306 individuals using standard laboratory procedures. Two different sequencing technologies were employed: Illumina HiSeqX Five/Ten was used to sequence 155 samples (conducted by Macrogen Japan Corporation [https://macrogen-japan.co.jp] in 2017 and 2018), and Illumina NovaSeq was used to sequence 156 individuals (by RIKEN sequencing platform in 2020). Both Illumina HiSeqX Five/Ten and NovaSeq platforms generated 150 bp paired-end reads and a mean depth of 18.4x. Out of the 311 samples sequenced, five were duplicates, sequenced with both platforms. Variant calling was performed using Dragen Bio IT-platform v3.5.7 (Illumina) with the GRCh38 human reference genome. Joint calling was performed with Dragen Bio IT-platform v3.6.3 (Illumina) and excluded low quality variants (QUAL<10.41 for SNPs and QUAL<7.83 for InDels) following Dragen default parameters. The joint called VCF files were then “lifted down” from the GRCh38 to the GRCh37 reference genome for downstream analyses.

In an attempt to understand the lack of correlation at the genus level between Ishida et al. (2020) study and ours, we utilized two methods to process fecal samples from five healthy subjects, in duplicate, for microbial analysis using 16S rRNA sequencing.

The first method, replicated from Ishida et al., involved preserving the fresh fecal sample in GuSCN solution, vortexing with glass beads, and treating with buffer-saturated phenol. After centrifugation, the supernatant was further extracted with phenol-chloroform and precipitated with isopropanol. The DNA, extracted by this bead-beating method, was then subjected to 16S rRNA sequencing. The resulting sequences were classified into operational taxonomic units (OTUs) using QIIME v2 (Bolyen et al., 2019) for analysis and Greengenes (DeSantis et al., 2006) for classification, with a 97% identity threshold.

The second method, used in our study, involved treating the fecal samples with lysozyme, achromopeptidase, and proteinase K for lysis, followed by phenol-chloroform separation and ethanol precipitation. The DNA, extracted by this enzymatic lysis method, was preserved at −80°C prior to 16S rRNA sequencing. The resulting sequences were grouped into OTUs using UCLUST, also with a 97% identity threshold.

We extracted rs671 variant (MAF 27%) from 275 individuals in our dataset. These individuals were categorized into two groups based on rs671 genotypes: major allele homozygous group (GG); and minor allele homozygous and heterozygous group (AG + AA). Additionally, we randomly selected 999 SNPs with allele frequencies ranging from 25% to 30% from our dataset and categorized them into two groups: major group (major allele homozygous group); and minor group (minor allele homozygous & heterozygous group), following the same approach. The fecal metabolite data was obtained from a previous study (Takeuchi et al., 2021; Takeuchi et al., 2023), and normalized using Blom normalization method, whereas cytokine and clinical data such as Homeostatic Model Assessment of Insulin Resistance (HOMA-IR) and BMI were normalized with inverse rank-based method.

Given the evidence that proinflammatory cytokines play a key role in modulating insulin signaling, our previous research identified 29 triangular relationships in which certain cytokines significantly affected IR markers, including HOMA-IR (Takeuchi et al., 2023). To investigate whether rs671 further shapes these cytokine–IR associations, we conducted causal mediation analyses for each of the 29 relationships using 1,000 SNPs (including rs671) employing R mediation package. For each SNP, we performed 2,000 mediation analyses (1,000 SNPs × 2 groups) per relationship (2,000 × 29 = 58,000 in total) across major and minor genotype groups, then calculated Z-scores based on the p-values for the Average Causal Mediation Effect (ACME), Average Direct Effect (ADE), and Total Effect (TE). With ACME measuring the portion of the total effect of a SNP on insulin resistance, that is mediated through cytokine levels, thus quantifying the indirect effect. ADE representing the portion of the total effect that is not mediated, reflecting the direct effect on insulin resistance, and TE the sum of ACME and ADE, representing the overall effect, as schematically represented in Supplementary Figure 2 . Because of occasional zero p-values, we added 0.001 before Z-transformation. We assessed the impact of rs671 by subtracting the Z-score of the minor group from that of the major group and evaluating significance.

Lastly, the rs671 association results were extracted from our GWAS analysis. Keeping only nominal associations between rs671 and microbial features, with beta values indicating the direction and magnitude of the associations.

To investigate the relationship between the host and gut microbiome in the Japanese population, we enrolled 306 individuals. Whole-genome sequencing (WGS) and fecal shotgun metagenome sequencing analyses were performed on these individuals (Materials and methods). After conducting per-sample and per-marker quality control for WGS, we obtained 296 samples and approximately 13 million variants ( Supplementary Figure 3 ). Additionally, we applied filtering based on prevalence and abundance, as well as transformation of the taxonomic compositional data (Materials and methods). This resulted in the identification of three phyla, six classes, seven orders, 10 families, 17 genera, 86 species and 98 bacterial pathways from the shotgun metagenome sequencing data ( Supplementary Table 1 ).

Despite the filtering and inverse normal transformation (Materials and methods), several taxa, particularly at the species level, did not exhibit a normal distribution ( Supplementary File 1 ). Instead, these taxa displayed a binary-like distribution, likely due to their low abundance in only a subset of individuals, which reflects the compositional nature of the microbiome at different taxonomic levels. Thus, by applying the Shapiro–Wilk test to the 227 taxa and bacterial pathways, we revealed that 140 taxa had W-statistics ≥ 0.95, indicating that their distributions approximated normality ( Supplementary Table 2 ). These taxa were thus suitable for analysis as quantitative phenotypes using linear regression models in GWAS. Conversely, 87 taxa had W-statistics below 0.95, failing to meet the normality criterion. These taxa displayed distribution patterns consistent with binary traits, warranting additional analysis through logistic regression models. The histograms of taxa distributions corroborated the Shapiro–Wilk test results, as taxa failing the normality threshold often exhibited skewed or bimodal distributions. This distinction ensured that each taxon was analyzed using the most appropriate statistical model, enhancing the reliability of the GWAS findings.

Next, we conducted GWAS between the QCed host genotypes and each filtered and transformed taxa and bacterial pathway (Materials and methods). GWAS results ( Table 1 ) revealed no significant associations at the phylum, class, order and species levels. However, we identified one significant (p < 5×10⁻⁸, genome-wide significance threshold) loci at the family level ( Figure 1A ), one at the genus level ( Figure 1B ), and eight significant loci at the bacterial pathway level ( Supplementary Figure 4 ). Nonetheless, when applying a more stringent threshold based on the number of tested phenotypes (p < 2.75×10^-11), all loci lost its significance.

The logistic regression GWAS conducted on the 87 binary taxa and bacterial pathways identified one intronic variant, rs4880022-C, located in the PAX5 gene, associated with the presence of the species unknown Clostridiales [meta_mOTU_v2_5805] at a p-value of 3.25×10⁻⁸ ( Table 1 , Figure 2 ). This association suggests that individuals carrying rs4880022-C have a higher likelihood of harboring this specific Clostridiales species (OR = 4.9, SE = 0.29). Additionally, several suggestive associations (p ≤ 1.0×10⁻⁵) were identified ( Supplementary Table 3 ). For example, a variant in unknown Clostridiales [meta_mOTU_v2_6852] on chromosome 19 (19:12646825:C:CA) showed a suggestive association with an OR of 0.236 (SE = 0.27, p = 9.79×10⁻⁸), while another variant in Clostridium innocuum [ref_mOTU_v2_0643] on chromosome 3 (3:16735689:G:T) was associated with an OR of 5.27 (SE = 0.32, p = 2.44×10⁻⁷). These suggestive associations may represent potential genetic influences on the presence of other low-abundance taxa, specially from the Clostridiales family, and warrant further investigation. However, after applying the stringent correction for multiple testing (p ≤ 9.32×10⁻¹¹), the association between rs4880022 and unknown Clostridiales [meta_mOTU_v2_5805] did not remain statistically significant.

We compared our results with the latest four studies on host–microbiome GWAS in Japanese (Ishida et al., 2020), European (Lopera-Maya et al., 2022; Qin et al., 2022), and multi-ancestry populations (Kurilshikov et al., 2021) ( Table 2 ). In our study, we analyzed 227 taxa and pathways from 296 individuals and found 11 independent significant loci. However, as observed in previous comparable studies, there was a lack of reproducibility in significant GWAS associations across different studies ( Supplementary Table 4 ), with a direction of effect consistent in only 46% of the associations, indicating that the observed associations may be due to random chance rather than true genetic effects. When considering the association analysis conducted in Japanese ethnicity by Ishida et al (Ishida et al., 2020), there was a lack of replicated significant associations between studies, with only four out of 38 associations being nominally significant in our results ( Supplementary Table 4 ). Additionally, the direction of effect for the significant associations between studies was inconsistent when compared to our summary statistics, with only 53% showing the same direction of effect. Overall, this inconsistency was demonstrated by a weak positive linear Pearson correlation and a p-value indicating weak correlation of 0.334 (p = 0.020) when comparing with beta values from all significant variants from previous studies. Moreover, a Pearson correlation of 0.374 (p = 0.021) was observed when comparing against Ishida et al. (2020) effect sizes from their suggested associations ( Supplementary Figure 5 ).

To further investigate the lack of reproducibility between studies, particularly those conducted on the same ethnicity (Japanese), we compared the raw core genus relative abundance composition from the study by Ishida et al. with the same genera’s relative abundance obtained from our 16S rRNA sequencing (Materials and methods). As shown in Figure 3 , the comparison revealed a low correlation (R = 0.38, p = 0.19) and a Bray–Curtis dissimilarity value of 0.46, implying a moderate dissimilarity between the two cohorts. This suggests that while the two cohorts share some common genera, they differ in their relative abundance composition.

In an attempt to understand the possible causes of these differences in relative abundances, we replicated the experimental workflow used in both studies to generate the 16S rRNA sequencing relative abundance data ( Table 3 ). We applied the same conditions and steps from both methods to five healthy Japanese samples (Materials and methods) and generated the relative abundance composition for each sample at the phylum and genus levels. First, we compared the replicates’ relative abundance at the genus level and found an average correlation between replicates of 97.4% when using the Ishida et al. method and 99.4% when performing our method ( Supplementary Figure 6 ). Then, when comparing between methods, our results showed an average correlation of 95.5% (p = 2.2 × 10⁻⁴) at the phylum level. Figure 4 illustrates a more straightforward comparison of relative abundance data, by extracting only the core genera from the method used by Ishida et al. and compared it with the relative abundance output from the method used in our study. Similar to our initial comparison of raw genus relative abundance composition depicted in Figure 3 , our analysis showed a higher relative abundance of Bacteroides and a decreased relative abundance of Bifidobacterium in the output from the Ishida et al. method. With an overall average correlation of core genus between the two methods of 77% (p = 4.12 ×10⁻⁶).

To overcome our limited sample size and find further associations between common host functional genetic variants and our taxonomic compositional data, we conducted a PheWAS across all taxonomic levels. We analyzed 1,412 high-impact variants, which likely cause a disruption in gene function, including nonsense mutations, frameshift mutations, and splice site mutations, with a minor MAF ≥ 0.05 (Methods, Supplementary Table 5 ). This approach allowed us to assess the impact of common, functional genetic variants on the relative abundance of gut microbial taxa.

Our PheWAS identified one significant frameshift variant located in the OR6C1 gene (rs5798345-CA, c.24dup, p.Glu9ArgfsTer10) associated with the relative abundance of Bacteroides uniformis (beta = 0.394, p ≤ 4.78 × 10⁻⁸; Table 4 ). The positive beta value indicates that individuals carrying the frameshift variant have a higher relative abundance of B. uniformis in their gut microbiome. The same frameshift variant demonstrated nominal significance in other Bacteroides species (Bacteroides fragilis, Bacteroides dorei, Bacteroidales_sp) and their corresponding taxonomic levels, all showing a consistent positive direction of effect ( Figure 5A ). Furthermore, we excluded Bacteroides uniformis and included other available Bacteroides species from our cohort that had been filtered out due to low prevalence and abundance. We then re-ran the PheWAS analysis using only rs5798345 to determine whether the higher taxonomic levels of B. uniformis retained their nominal significance, as well as to identify any additional Bacteroides species with nominal significance and consistent direction of effects. Notably, the higher taxonomic levels retained their nominal significance, and we identified four additional nominally significant independent species (B. stercoris, B. thetaiotaomicron, B. massiliensis, and B. plebeius) with concordant directions of effect ( Figure 5B ; Supplementary Table 6 ). These findings suggest that the OR6C1 variant may have a broader impact on the abundance of taxa within the Bacteroides genus.

To enhance the identification of loci, specifically genes, that may be linked to gut microbiome composition in our cohort, we employed the strategy of mapping our high-impact-annotated variants to their corresponding genes. Subsequently, we conducted a gene-based analysis using MAGMA software (Materials and methods). However, we were unable to find any robust functional associations (p < 2.29×10⁻⁶; Supplementary Table 7 ).

According to previous studies, causal mediation analyses demonstrated that some inflammatory cytokines may mediate the effects of fecal carbohydrates on insulin resistance (Takeuchi et al., 2023). Here, we investigated whether host genetics—particularly rs671—could further shape these relationships in our cohort. This SNP, which profoundly affects alcohol metabolism, has been linked to cardiovascular diseases, cancer (Chang et al., 2017; Zhang et al., 2023; Koyanagi et al., 2024), and susceptibility to type 2 diabetes in males (Spracklen et al., 2020).

Our GWAS results indicate that the rs671 variant is nominally associated with the bacterial chemotaxis pathway [PATH:ko02030], exhibiting a positive beta value ( Supplementary Table 8 ). Additionally, rs671 shows nominal associations with the relative abundance of Clostridium innocuum and Streptococcus salivarius, as well as their higher taxonomic levels. For S. salivarius, associations extend to the genus Streptococcus, family Streptococcaceae, order Lactobacillales, and class Bacilli. Both species are part of the phylum Firmicutes. The consistent positive beta values across these taxa suggest a potential relationship between the rs671 minor allele and increased abundance of these bacteria, as well as enhanced representation of bacterial chemotaxis functions within the gut microbiota, though the statistical significance of these associations were lost after multiple testing correction.

To assess the impact of rs671 on the mediation of fecal carbohydrates’ effects on insulin resistance, we calculated Z-score differences for the ACME, ADE, and TE across 29 cytokine–IR relationships for 1,000 SNPs ( Supplementary Table 9 ; Materials and methods). Our results revealed one pair (ACME), six pairs (ADE), and one pair (TE) that surpassed a Z-score difference of 2, suggesting that this EAS-specific variant may influence how carbohydrates affect IR. These findings underscore the potential importance of rs671 in modulating gut microbial pathways and cytokine-mediated IR processes.