The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Ethics Committee of Mae Fah Luang University (EC23251-11) on 19 December 2023. The fecal samples analyzed in this study were obtained from a previous study (Gruneck et al., 2022), in which written informed consent was obtained from all participants. The current analysis used a subset of those previously collected and anonymized fecal samples, without involving any new participant contact or data collection.

This study included 67 school-aged children (6–15 years old), selected as a subset of 127 children from Ban Huai Rai Samakee Elementary School in Chiang Rai, Thailand (20 °16′13.6″N 99 °49′43.4″E). Details regarding participant recruitment and BMI classification have been previously described (Gruneck et al., 2022). Briefly, recruitment was conducted through voluntary participation coordinated by the school administration, with informed consent obtained from parents prior to participation. Inclusion criteria required children to be healthy, free from acute gastrointestinal diseases (e.g., gastroenteritis, diarrhea), and not to have taken antibiotics within 2 months before specimen collection. The number of participants was reduced from our previous study to balance group sizes and minimize confounding effects of ethnicity on gut microbiota composition by equalizing the number of individuals from each ethnic group (Gruneck et al., 2022). Participants were classified into normal weight (N = 30), overweight (OW = 20), and obese (OB = 17) groups based on gender-specific BMI-for-age z-scores, using World Health Organization cut-offs for children aged 5–19 years (De Onis et al., 2007). Fecal samples were collected on site in sterile containers, stored in ice boxes, and immediately transported to the laboratory for storage at −80 °C.

Fecal microbial DNA was extracted following the same protocol described in our previous study (Gruneck et al., 2022). For library preparation and NGS analysis, the extracted genomic DNA samples were submitted to Novogene Co., Ltd. (Singapore) and amplified using primers 341F-5′ CCTAYGGGRBGCASCAG 3′ and 806R-5′ GGACTACNNGGGTATCTAAT 3′, targeting V3-V4 variable regions of 16 s rRNA gene. Pooled DNA library was quantitatively assessed using Qubit® 2.0 (Life Technologies, CA, United States), before loading into Illumina Novaseq Xplus (Illumina, San Diego, CA, United States) paired-end reads 250 × 2 bp according to the standard protocols.

Raw sequencing data was processed using Quantitative Insights Into Microbial Ecology 2 (QIIME2) pipeline version 2024.5.0 (Bolyen et al., 2019). Briefly, quality control was performed using the Cutadapt plugin (Martin, 2011) to remove adapter and primer sequences from both sequencing reads. Quality filtering, denoising, and merging were conducted using Divisive Amplicon Denoising Algorithm2 (DADA2) plugin (Callahan et al., 2016), generating amplicon sequence variants (ASVs) for phylogenetic and taxonomic analyses. Multiple sequence alignment of representative ASV sequences was performed using the MAFFT algorithm (Katoh and Standley, 2013). Phylogenetic tree was constructed using a maximum-likelihood method implemented in FastTree version 2.1.3 (Price et al., 2009). Taxonomic classification of ASVs was performed using the classify-sklearn Naive Bayes taxonomic classifier (Bokulich et al., 2018), with taxonomic assignment based on the Silva ribosomal RNA gene database (silva-138.2-ssu-nr99, released July 11, 2024).

Downstream analyses were conducted in R software (version 4.4.3) (R Core Team, 2025), with detailed summaries provided in Supplementary Data S1. Alpha and beta diversity were assessed at a rarefaction depth of 48,426 reads per sample, encompassing a total of 1,563 taxa (ASVs) (Supplementary File S1), using the phyloseq package (version 1.50.0) (McMurdie and Holmes, 2013). Alpha diversity metrics included Observed Species, Chao1 (species richness), Shannon (species richness and evenness), and PD Whole Tree (phylogenetic diversity). Group comparisons were performed using the Kruskal–Wallis test and visualized using ggplot2 (version 3.5.2) (Wickham, 2016). Beta diversity was calculated using Unweighted UniFrac, Weighted UniFrac, and Bray–Curtis distance metrics. Principal Coordinate Analysis (PCoA) was performed using the ordinate function (phyloseq) and visualized with plot_ordination (based on ggplot2). Permutational Multivariate Analysis of Variance (PERMANOVA) based on the three distance metrics was conducted using the adonis2 function in the vegan package (version 2.6–10) (Oksanen et al., 2025), with 999 permutations. Pairwise PERMANOVA (permutation = 999) was then performed using the pairwise.adonis function in the pairwiseAdonis package (version 0.4.1). Homogeneity of group dispersions (permutation = 999) was tested using the betadisper function in vegan. Gut microbiota abundance was compared using the Kruskal–Wallis test (p < 0.05), followed by Dunn’s test with Benjamini–Hochberg (BH) correction (hereafter referred to as the q-value). The relative abundance of gut microbiota at the genus level was visualized using the ComplexHeatmap package (version 2.22.0) (Gu et al., 2016; Gu, 2022). Hierarchical clustering of rows and columns was performed based on Spearman correlation.

Differentially abundant gut microbiota at the genus level (26 variables, CLR-transformed), positive-ion metabolites (95 variables, log₂-transformed), and negative-ion metabolites (38 variables, log₂-transformed) across BMI groups were integrated using Multiple Factor Analysis (MFA) to identify key features contributing to overall variation. MFA was conducted using the FactoMineR package (version 2.11) (Lê et al., 2008), and the factoextra package (version 1.0.7) (Kassambara and Mundt, 2020) was used to visualize variable contributions and correlations in the MFA space. To further investigate the influence of BMI groups on variation in gut microbiota (144 genera) and metabolite profiles (positive-ion mode = 95; negative-ion mode = 38), as well as the relationships between them, redundancy analysis (RDA) was performed using the vegan package. The BMI group was treated as a constrained explanatory variable, while the gut microbiota (CLR-transformed) and metabolite data (log₂-transformed and scaled) were used as response variables. The significance of the constraints was assessed using an ANOVA-like permutation test (999 permutations). RDA biplots were generated with type 1 scaling (distances between samples represent similarity) and type 2 scaling (angles between arrows reflect correlations among variables), and visualized using ggplot2.

The functional profiles of gut microbiota were predicted using PICRUSt2 (version 2.6.2) based on KEGG Orthology (KO). Annotated genes were represented as K numbers and mapped against the KEGG database (http://rest.kegg.jp/get/br:ko00001, accessed May 30, 2025) to obtain pathway-level functional annotations. Differentially abundant KOs across BMI groups were identified using the ALDEx2 package (version 1.38.0). Welch’s t-test (we.ep) and Welch’s t-test Benjamini-Hochberg adjusted q-values (we.eBH) were used to determine significance. Pairwise comparisons included N vs. OW, N vs. OB, and OB vs. OW. Associations were evaluated between gut microbiota (144 genera) and 5,819 predicted KOs, between fecal metabolites (positive-ion mode: 95; negative-ion mode: 38) and 5,819 predicted KOs, and between fecal metabolites and KO abundance stratified by contributing ASVs at the genus level (KO|Genus) based on significantly differentially abundant KOs identified via ALDEx2. All associations were assessed using the Hierarchical All-against-All Association (HAllA) method (Ghazi et al., 2022), implemented in Python (version 3.10.16). Significant associations (Spearman correlation coefficient r ≥ |0.7|, q < 0.05) were further analyzed through network analysis using the igraph package (version 2.1.4) (Csárdi et al., 2025). Community detection in the network was performed using the Louvain clustering algorithm.

Microbial diversity was assessed at a sequence depth of 48,426 reads/sample, retaining 1,563 ASVs. Alpha diversity analysis showed no significant difference among BMI groups (Supplementary Figure S1) based on four indices including observed species, Chao1, Shannon, and PD. Beta diversity based on the unweighted UniFrac (F_{paired. PERM} = 1.95, q = 0.006) revealed significant differences in community membership between N and OB groups, with no significant difference in dispersion (Supplementary Figure S2). No significant differences in overall community structure and composition were observed using the weighted UniFrac and Bray-Curtis distance matrices.

Based on 1,563 ASVs conserved across all samples, 11 known ASVs were identified at the phylum level, 17 at class, 30 at order, 48 at family, and 143 at genus level. ASVs that could not be classified from class to genus were grouped as unclassified. Microbiome composition analysis revealed that the five most abundant phyla across all BMI groups were Bacteroidota, Bacillota, Actinomycetota, Verrucomicrobiota, and Fusobacteriota (Supplementary Figure S3A). A significant difference was observed for Actinomycetota (q = 0.02), Verrucomicrobiota (q = 0.02) and Elusimicrobiota (q = 0.03) between N and OB groups; however, the latter phylum was present in only one sample within the OB group (Supplementary Figure S3B; Supplementary File S3).

At the class level, three bacterial taxa showed significant differences between BMI groups (Supplementary Figure S4A). Coriobacteriia and Negativicutes were significantly more abundant in the OB group compared to both N and OW groups (q < 0.05), while Elusimicrobia was also more abundant in OB compared to N (q = 0.023). Other dominant classes such as Bacteroidia, Clostridia, Verrucomicrobiia, Actinobacteria, Bacilli, Fusobacteriia, and Alphaproteobacteria did not show significant differences among BMI groups.

At the order level, Coriobacteriales, Veillonellales-Selenomonadales, Elusimicrobiales, Acidaminococcales, and Caulobacterales showed the highest abundance in the OB group compared to the N and OW groups (q < 0.05) (Supplementary Figure S4B). At the family level, Rikenellaceae had the lowest abundance in the OB group compared to N and OW (q < 0.05) (Supplementary Figure S5). Barnesiellaceae and Marinifilaceae were more abundant in the N group than in OB (q < 0.05), whereas Coriobacteriaceae, Selenomonadaceae, Acidaminococcaceae, and Caulobacteraceae were most enriched in the OB group relative to the N and OW groups (q < 0.05). Elusimicrobiaceae abundance was significantly higher in OB than in N (q = 0.008).

At the genus level, Segatella and Bacteroides were the most prevalent genera across all groups (Supplementary Figure S6). Among the 24 genera with significant differences between BMI groups, five were most prevalent in the OB groups (q < 0.05), including Faecalibacterium, Collinsella, Megamonas, Brevundimonas, and Phascolarctobacterium (Figure 1). Additionally, Massiliprevotella, Tractidigestivibacter, and Elusimicrobium were also more abundant in OB compared to N (q < 0.05). Prevotellaceae UCG-003 was most enriched in the OW group, while the abundance of Howardella was also higher in OW compared to N (q < 0.05). In contrast, Lachnospiraceae UCG-010, Eisenbergiella, Ruthenibacterium, and Paludicola exhibited lowest abundance in the OB group (q < 0.05). The abundance of 10 genera, including Alistipes, Barnesiella, Hungatella, Butyricimonas, Odoribacter, Intestinimonas, Anaerotruncus, Holdemania, Dielma, and GCA-900066755, were also significantly less abundant in OB compared to N (q < 0.05). Overall, the gut microbiota composition in OB participants was more distinct from that of N than OW, highlighting obesity-associated changes in microbial profiles.

Using a PLS-DA cutoff of VIP > 1, a total of 133 differential fecal metabolites were identified across BMI groups (Figures 2A,B), of which 95 were detected in positive-ion mode and 38 in negative-ion mode. The mean log₂ fold change (log₂FC) of positive-ion metabolites relative to the N group revealed that 22 positive-ion features were upregulated, 70 were downregulated, and 3 remained unchanged in the OW group (Supplementary Figure S7; Supplementary File S4). In the obese (OB) group, 40 positive-ion features were upregulated, while 55 were downregulated. For negative-ion metabolites, 31 features were consistently upregulated in both OW and OB groups, whereas only 7 features were downregulated in these two groups compared to the N group. Additionally, direct comparisons of log₂-transformed positive-ion metabolite abundances revealed that 27 features were most enriched in the OB group (q < 0.05), whereas 30 and 4 features were most abundant in the N and OW groups, respectively (q < 0.05; Supplementary File S4). For negative-ion metabolites, 30 out of 38 features were least abundant in the N group compared to both OW and OB groups (q < 0.05), suggesting an increased abundance of these metabolites in individuals with higher BMI. Notably, the OW group exhibited a negative metabolite profile more similar to that of the OB group than to the N group. Nearly half of these negative-ion metabolites (16 features) were significantly enriched in the OB group (q < 0.05), further emphasizing their association with higher BMI.

Further sPLS-DA analysis based on a tuned model (VIP score > 1) identified 70 and 23 metabolite features (across both ion modes) as important for discriminating between BMI groups in components 1 and 2, respectively (Supplementary File S4). Component 1, which explained 53% of the variance, revealed fewer positive-ion metabolites associated with the OB group (21, AUC = 96, p < 0.0001) compared to the N group (23, AUC = 97, p < 0.0001), while 22 negative-ion metabolites contributed to the separation of OB from N (Figure 2C). The top five metabolites enriched in OB were neg4 (Ala Lys Pro Gln), pos53 (unclassified), neg11 (JWH 018 N-(5-hydroxypentyl) metabolite-d5), neg8 (Sarcostin), and neg 21 (15-keto-Prostaglandin E2). Only four negative metabolites contributed to the N group. In component 2, which accounted for for 14% of explained variance, more positive-ion compounds were associated with OB (16, AUC = 0.94, p < 0.0001) than with OW (7, AUC = 1, p < 0.0001) (Supplementary Figure S8B). The top contributors for OB included pos20 (L-Phenylalanyl-L-valyl-L-leucyl-L-prolyl-L-tryptophyl-N ~ 5 ~ −(diaminomethylidene)-L-ornithyl-L-isoleucine), pos88 (Amaranthussaponin II), pos65 (PC(16:0/9:0(CHO))), pos48 (unclassified), and pos23 (unclassified), while pos37 (unclassified) and pos91 [Endothall, di(N, N-dimethyltridecylamine) salt (1:2)] were highly associated with OW. Among these, pos88 (Amaranthussaponin II; C₄₈H₇₄O₂₀; HMDB0041352) is an organic compound classified as a lipid-like molecules and functions as an energy source; however, no specific KEGG pathway annotation was identified.

Multiple factor analysis (MFA) integrating gut microbiota abundance and fecal metabolite profiles revealed that inter-individual variations across BMI groups were primarily explained by metabolite profiles on Dimension 1 (34.2% explained variance), and by gut microbiota composition on Dimension 2 (10.5% explained variance) (Figure 3; Supplementary File S5). On Dimension 1, both positive-ion (cos² = 0.76) and negative-ion (cos² = 0.82) metabolites were major contributors to the observed variations. This was best illustrated by a distinct separation between the OB group (coordinate = 1.60, p < 0.0001) and the N group (coordinate = −1.77, p < 0.0001), reflecting contrasting metabolite profiles, especially in negative-ion mode (Supplementary Table S1). Several negative-ion compounds, including neg8, neg22, neg21, neg11, neg4, neg6, and neg19, were highly correlated with this dimension and enriched in the OB group. In contrast, neg17, neg18, neg30, and neg24 were associated with the N group. Although gut microbiota showed a high contribution (79.74%) to dimension 2, no single dimension (Dim 1 to Dim 4) adequately represented microbiota composition (summed cos² < 0.5). Overall, the MFA demonstrated that BMI was more strongly associated with variations in fecal metabolite profiles than with gut microbiota composition in this study.

Furthermore, RDA revealed that BMI groups (X) significantly influenced metabolite profiles (Y), with the model yielding a high canonical correlation (F_anova.cca = 30.33, p = 0.001) and a constrained proportion of 53.84% (variance in Y explained by X; adjusted R² = 0.52). When gut microbiota composition at the genus level was included in the RDA model of metabolites, the model remained significant (F_anova.cca = 4.92, p = 0.001), though the constrained proportion dropped substantially to 15.93% (adjusted R² = 0.13). In the combined model, BMI explained 75.49% of the variation in RDA1 (F_anova.cca = 7.44, p = 0.001) and 24.51% in RDA2 (F_anova.cca = 2.41, p = 0.001). Distinct clustering of BMI groups was observed on scaling 1 (Figure 4A), where both metabolite features and gut microbiota showed an alignment with BMI group centroids, as indicated by arrow directionality. Among the top 100 contributing features, several positive-ion metabolites and bacterial genera, including Collinsella, Romboutsia, and Meganomas, shared vector directions in the OB group space and clustered near the group centroid. Conversely, some positive-ion metabolites and gut microbiota such as Paraprevotella, Akkermansia, and Odoribacter were clustered near the N group centroids. Many of the negative-ion metabolites also exhibited vector directions aligned with the OB and OW group spaces, while positive-ion features were enriched in the N group space. On scaling 2 (Figure 4B), both Romboutsia and Collinsella, associated with the OB group, shared directional similarity and showed strong correlations with several positive-ion metabolites including pos87, pos75, pos76, and pos69 (Supplementary Table S2). In contrast, Megamonas demonstrated strong correlations with pos68, pos67, pos64, pos88, pos65, pos48, pos49 in the OB group space. Within the N group space, Paraprevotella showed a strong correlation exclusively with positive-ion metabolites, whereas Akkermansia, positioned further from the N group centroid, correlated strongly with neg24 (Asn His Val), neg30 (3′?-Isopravastatin), Odoribacter, neg17 (1-[2-Fluoro-4-(methanesulfonyl)phenyl]-1,4-diazepane), and neg18 (Hoechst 33342). Collectively, these RDA results highlight a strong influence of BMI on fecal metabolite profiles, as well as associations between metabolite features and specific gut microbiota that varied by BMI group. Moreover, the relatively lower explanatory power of gut microbiota (genus level) in the integrated model indicated fewer interactions between these two dataset in this study.

We further examined the associations between metabolites (positive- and negative-ion modes) and gut microbiota at the genus level using HAllA. The analysis revealed no significant association patterns between these two datasets within each BMI group after p-value adjustment.

Gut microbiota functional prediction identified a total of 5,819 KEGG genes across 54 KEGG level 2 pathways and 420 KEGG level 3 pathways (Supplementary File S6). Among these, K03088 (rpoE; RNA polymerase sigma-70 factor, ECF subfamily) was the most abundant KO across all BMI groups (Supplementary Figure S9). The top five most abundant level 2 pathways included protein families: signaling and cellular processes (16.77%), protein families: genetic information processing (10.06%), carbohydrate metabolism (8.06%), protein families: metabolism (6.93%), and amino acid metabolism (5.62%). At KEGG level 3, the most abundant pathways were transporters (9.16%), two-component system (4.75%), enzymes with EC numbers (4.14%), ABC transporters (3.53%), and function unknown (2.64%). Differential abundance analysis of KEGG orthologs (KOs) between BMI groups using ALDEx2 revealed 18, 9, and 17 differentially abundant KOs for N vs. OW, N vs. OB, and OB vs. OW, respectively (Supplementary Figures S10, S11). However, only differentially abundant KOs identified in N vs. OB comparisons remained statistically significant after p-value adjustment (we.eBH < 0.05). Mapping these differentially abundant KOs (we.ep < 0.05 or we.eBH < 0.05) to KEGG pathways revealed distinct enrichment patterns across ALDEx2 comparisons. In the N vs. OW comparison, protein families: signaling and cellular processes (level 2) and transporters (level 3) were the most enriched pathways (Figure 5A). In the N vs. OB comparison, protein families: signaling and cellular processes (level 2) was the most enriched (Figure 5B), while in the OB vs. OW comparison, the top three enriched pathways were oxidative phosphorylation (level 3), energy metabolism (level 2), and carbohydrate metabolism (level 2) (Figure 5C).

We further explored the associations between gut microbiota at the genus level and predicted functions (KOs) using HAllA. The analysis revealed dynamic shifts in association patterns across BMI groups (Supplementary Figure S12; Supplementary File S7). As BMI increased (N ⇢ OW ⇢ OB), the number of significant association clusters between gut microbiota and KOs declined markedly (N = 8,767; OW = 2,281; OB = 1,934). For example, the number of KOs significantly associated with major genera (q < 0.05, regardless of correlation strength), such as Segatella and Bacteroides, substantially decreased with increasing BMI (N = 1,709 and 927; OW = 523 and 523; OB = 314 and 120, respectively). Segatella also maintained strong associations with several abundant KOs ranked among the top 20 (Supplementary Figure S13). In contrast, OB-enriched genera including Faecalibacterium (OB = 191; N = 29), Megamonas (OB = 131; N = 8), Collinsella (OB = 74; N = 41), and Brevundimonas (OB = 59; N = 0) exhibited a notable increase in the number of associated KOs based on strong and significant HAllA clusters (r ≥ |0.7|, q < 0.05; Supplementary Figure S14). Conversely, Phascolarctobacterium showed relatively few KO associations in the OB group (OB = 5) compared to the N group (N = 245). Our results suggested a major reorganization in the functional interactions of the gut microbiome in obesity, characterized by fewer significant associations among dominant genera and increased coordination between OB-enriched genera and KOs. To visualize unique and shared associations across groups, we extracted significant and strong correlation clusters (excluding unclassified genera) from HAllA outputs (Spearman correlation coefficient r ≥ |0.7|, q < 0.05). This yielded 2,837 unique associations in the N group, 4,360 in OW and 4,234 in OB. Fewer associations were shared between the N and OB groups (1,160) compared to the N and OW groups (1,987), while 1,436 associations were shared between OW and OB. Among the top unique associations, Segatella was the dominant genus in the N group, showing strong negative correlations with most KOs. In contrast, Ruthenibacterium in the OW group demonstrated predominantly strong positive associations with KOs. In the OB group, genera such as Blautia, Segatella, NK4A214 group, and Faecalibacterium exhibited exclusively positive correlations with KOs. Additionally, OB-enriched genera like Collinsella and Brevundimonas also showed solely positive associations, whereas Megamonas displayed predominantly negative correlations with KOs. We also found that only one significantly different KO identified by ALDEx2 in the N vs. OB comparison overlapped with genus–KO associations detected by HAllA. This overlapping KO was K08281 (pncA; nicotinamidase/pyrazinamidase [EC:3.5.1.19, 3.5.1.-]), which was less abundant in the OB group. In contrast, greater overlap was observed in the N vs. OW (5 KOs) and OW vs. OB (4 KOs) comparisons. These findings suggest a progressive reduction in the connection between differentially abundant KOs and their associations with gut microbiota as BMI increases. In other words, changes in functional profiles may not directly correspond to specific microbiota–function relationships, particularly in higher BMI groups.

To further investigate the interactions between gut microbiota and predicted functions (KOs), we performed network analysis using the strong HAllA correlation clusters for each BMI group. Louvain clustering was applied to detect community structures and connection patterns between gut microbiota and KOs. The analysis revealed increased network modularity and decreased centralization in the OB group compared to the OW and N groups, indicating more distinct community structures than centralized interactions within the OB network (Supplementary File S7). These findings implied that in N and OW individuals, microbial functions were more diffusely interconnected, whereas in obese individuals, functional interactions became more compartmentalized and structured. Segatella and Bosea served as a dominant hub in the N and OW groups, respectively, each displaying the highest connectivity to KOs. However, in the obese (OB) group, the hub centrality of both Segatella and Bosea declined by more than half, and they were surpassed by Megamonas, which became the new dominant hub. This marked shift suggested a substantial reorganization of the gut microbiota–function interaction network in obesity, with central functional roles redistributed among different taxa.

Determining the associations between metabolites and predicted functions (KOs) using HAllA revealed dynamic shifts in association patterns across BMI groups (Supplementary Figure S15). The number of significant association clusters (q < 0.05) between metabolites and KOs declined with increasing BMI, with 467 clusters in the N group, 13 in the OW group, and 181 in the OB group. When focusing on strong associations (r ≥ |0.7|, q < 0.05), 22 out of 54 metabolites were enriched by more than 200 KO associations in the N group. Among these, neg34 and pos84 were the top two metabolites with the highest number of KO associations, each linked to 322 associations (Supplementary Table S3). In the OW group, pos30 and pos43 were the dominant KO-associated metabolites (Supplementary Figure S16A). Notably, pos30 also remained predominant in the OB group (Supplementary Figure S16B), associated with 169 KOs. In addition, all metabolite–KO associations in the OW group were exclusively positive. Our results suggested a decline in functional interactions between metabolites and KOs in overweight individuals compared to the normal group.

Furthermore, none of the differentially abundant KOs identified by ALDEx2 overlapped with the strong metabolite–KO associations detected by HAllA. By integrating significantly strong HAllA associations, differentially abundant metabolites between BMI groups, and discriminant features from the sPLS-DA tuned model (component 1; VIP > 1) revealed that OB-discriminant metabolites showed a high number of associations with KOs in the N group, whereas N-discriminant metabolites predominantly interacted with KOs in the OW and OB groups (Figure 6). This pattern suggested a contrasting enrichment profile of metabolite–KO associations across BMI groups, particularly between the N and OB groups.

Integrating metabolites and the stratified output from PICRUSt2 (KO abundance stratified by contributing ASVs at the genus level, KO|Genus) into HAllA revealed that the number of significant association clusters (q < 0.05) between metabolites and KO|Genus substantially declined in the N group compared to the higher BMI groups, with 6 clusters in the normal (N) group, 154 in the OW group, and 91 in the OB group (Supplementary File S9). In the OW group, pos17 (PG(22:6(4Z,7Z,11E,13Z,15E,19Z)-2OH(10S,17)/i-20:0)) was the major metabolite forming the highest number of clusters with KO|Genus (109 clusters), while pos23 (unclassified) and pos93 (Fumitremorgin B) shared 17 clusters. In the OB group, pos2 (Riboflavin, Vitamin B2) was the dominant cluster-associated metabolite (57 clusters), while pos30 formed 6 clusters. Both pos17 and pos2 served as major hubs in the interaction networks of the OW and OB groups, respectively. More than 85% of the significant associations were strong across BMI groups. Further analysis of replicated associations based on shared KO|Genus patterns (i.e., KO|Genus groups that appeared with the same set of metabolites) revealed that only K00179|Collinsella, an OB-enriched genus, formed a strong positive association with pos30 (r = 1, q < 0.0001), and this association was unique to the OB group. None of the other OB-enriched genera, including Faecalibacterium, Brevundimonas, Megamonas, and Phascolarctobacterium, were KO-contributing taxa linked to any metabolites. Most genera (KO|Genus) associated with metabolites were not among the significantly abundant taxa across BMI groups. In addition, none of the strong associations were shared between the N and OW groups or between the N and OB groups. Twenty-eight associations were shared between the OW and OB groups, all involving pos30, pos43, and pos95 (Ala Phe Pro Glu). Specifically, pos30 and pos95 were associated with KO-contributing Bacteroides, while pos43 interacted with both KO-contributing Bacteroides and UCG-002. None of the differentially abundant KOs identified by ALDEx2 overlapped with the strong metabolite–KO|Genus associations detected by HAllA. Integration of significantly strong HAllA associations, differentially abundant metabolites between BMI groups, and discriminant features from the sPLS-DA tuned model (component 1; VIP > 1) yielded a similar pattern to the HAllA metabolite-KO associations. OB-discriminant metabolites were predominantly associated with KO|Genus in the N group, while N-discriminant metabolites showed more frequent associations with KO|Genus in the OW and OB groups (Supplementary Figure S17).

Collectively, our study demonstrated that gut microbiota interacted with metabolites primarily through functional bridging mechanisms rather than through direct associations. Interaction profiles involving metabolite variables were more similar between the higher BMI groups (OW and OB), highlighting a shared pattern of microbiome-metabolite crosstalk in individuals with elevated BMI. Significantly abundant features (gut microbiota or KOs) contributed less to the metabolite–KO or metabolite–KO|Genus interaction space across BMI groups. Pos30 served as a key metabolite in both OW and OB groups. Its strong association with the OB-enriched genus, Collinsella, via K00179 (iorA; indolepyruvate ferredoxin oxidoreductase, alpha subunit [EC:1.2.7.8]) appeared to represent an OB-specific functional signature in this study.