Gut microbiota dysbiosis and metabolic perturbations of bile/glyceric acids in major depressive disorder with IBS comorbidity

From October 2021 to August 2023 at the Mental Health Center of West China Hospital, Sichuan University, we enrolled 120 patients diagnosed with MDD using the Mini International Neuropsychiatric Interview (27), assessed by senior psychiatrists. In addition, we recruited 70 healthy subjects without MDD and IBS as the control group (healthy control [HC]). To be eligible for participation in this study, participants had to meet the following criteria: (i) age between 18 and 65 years (inclusive), (ii) education level of primary school or higher with a sufficient understanding of the study content, (iii) voluntary consent to participate in the study after being informed about its purpose and procedures, and finally, (iv) no medication use for more than 2 weeks or intermittent medication use for less than 3 days within a 2-week period. However, individuals were excluded if they presented with severe physical illness, organic brain disease, neurological disorder, or cognitive impairment. Moreover, individuals with comorbid mental illnesses such as schizophrenia and bipolar disorder were also excluded. Furthermore, individuals who consumed probiotics daily prior to their arrival at our center were also excluded from participating in this study (Fig. 1).

After arrival at our center, all subjects completed the following battery of tests: IBS Severity Scoring System (IBS-SSS) (28), Hamilton Depression Scale (HAMD-17) (29), Hamilton Anxiety Scale (HAMA-14) (30), Screening Depression Questionnaire (PHQ-9) (31), and Generalized Anxiety Disorder Scale (GAD-7) (32). The four scales used in this study (HAMA-14, GAD-7, PHQ-9, and HAMD-17) collectively constituted the participants' emotional state.

The IBS-SSS employs a questionnaire to assess the severity of abdominal pain, frequency of abdominal pain, degree of abdominal distension, defecation patterns, and impact on patients' quality of life. Scores are assigned using a visual scoring method ranging from 0 to 500, with higher scores indicating more severe symptoms. Based on test scores, participants were categorized into three levels: mild (75–175 points), moderate (176–300 points), and severe (301–500 points). MDD patients with an IBS-SSS score <75 were classified as the MDD without IBS group (MDD); conversely, those with an IBS score ≥75 were classified as the MDD with IBS group (MDD with IBS).

Demographic data were also collected for all subjects, including age, sex, smoking, drinking, body mass index (BMI), religion, civil status, employment status, and monthly household income. Meanwhile, stool and serum samples were also collected on the day of subject enrollment. Among them, 71 participants lacked a biological sample (Fig. 1).

DNA was extracted from fecal samples, and the purity and integrity of the extracted DNA were assessed by agarose gel electrophoresis. DNA concentration was accurately quantified using a Qubit fluorometer (Thermo Fisher Scientific). Subsequently, library preparation was performed, which included DNA fragmentation, end repair, A-tailing, adapter ligation, purification, and PCR amplification. After library construction, library quality and fragment size distribution were evaluated using an Agilent 2100 Bioanalyzer and quantitative PCR (qPCR). High-throughput sequencing was conducted on the Illumina NovaSeq 6000 platform using a paired-end 150 bp (PE150) mode.

Raw sequencing data underwent quality control with fastp (v0.19.3), including removal of adapter contamination, low-quality reads, and reads shorter than the threshold, resulting in high-quality clean reads. Clean reads were then aligned to the human reference genome GRCh38 (RefSeq: GCF_000001405.40↗) using Bowtie2 (v2.3.4) to filter out host-derived contamination. The remaining non-host reads were assembled both individually and in a combined manner using MEGAHIT (v1.2.9). Assembly quality was assessed with QUAST, considering metrics such as N50, total number of contigs, and total assembly length. Contigs of length ≥ 500 bp were subjected to open reading frame prediction using MetaGeneMark (v3.38). Predicted genes were filtered, standardized in naming, and clustered at 95% sequence similarity using CD-HIT (v4.8.1) to generate a non-redundant gene catalog (Unigenes). Clean reads from each sample were mapped back to the Unigene sequences using Bowtie2 to obtain read counts per gene. Unigenes with mapped reads ≤2 across all samples were removed to produce a high-quality gene set for downstream analysis. For taxonomic annotation, the protein sequences corresponding to Unigenes were aligned against bacterial, fungal, archaeal, and viral protein sequences extracted from the NCBI NR database (version 2022.05) using DIAMOND (v6.24.20). Taxonomic classification was performed based on the lowest common ancestor algorithm implemented in MEGAN (v0.9.24), generating species abundance profiles from kingdom to species levels. For functional annotation, DIAMOND was used to align Unigene protein sequences against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (version 2022.05). For each sequence, the annotation with the highest bit score >60 was retained for subsequent KEGG pathway annotation and functional analyses.

All of the standards of targeted metabolites were obtained from Sigma-Aldrich (St. Louis, MO, USA), Steraloids Inc. (Newport, RI, USA), and TRC Chemicals (Toronto, ON, Canada). All the standards were accurately weighed and prepared in water, methanol, sodium hydroxide solution, or hydrochloric acid solution to obtain individual stock solution at a concentration of 5.0 mg/mL. An appropriate amount of each stock solution was mixed to create stock calibration solutions. Formic acid was of Optima grade and obtained from Sigma-Aldrich (St. Louis, MO, USA). Methanol (Optima LC-MS), acetonitrile (Optima LC-MS), and isopropanol (Optima LC-MS) were purchased from Thermo-Fisher Scientific (Fair Lawn, NJ, USA). Ultrapure water was produced by a Mill-Q Reference system equipped with an LC-MS Pak filter (Millipore, Billerica, MA, USA).

Serum samples were thawed in an ice bath to minimize degradation. A 20 µL aliquot of plasma was transferred into a 96-well plate, followed by the addition of 120 µL of ice-cold methanol solution containing partial internal standards. After vigorous vortexing for 5 min, the mixture was centrifuged at 4,000 × g for 30 min. The plate was then transferred to an Eppendorf epMotion Workstation (Eppendorf Inc., Hamburg, Germany). Subsequently, 30 µL of the supernatant was transferred to a new 96-well plate, and 20 µL of freshly prepared derivatization reagent (200 mM 3-NPH in 75% aqueous methanol and 96 mM EDC-6% pyridine solution in methanol) was added to each well. The plate was sealed and incubated at 30°C for 60 min to allow derivatization. After the reaction, 330 µL of ice-cold 50% methanol solution was added to each well for dilution. The plate was then placed at −20°C for 20 min, followed by centrifugation at 4,000 × g for 30 min at 4°C. Subsequently, 135 µL of the supernatant was transferred to a new 96-well plate preloaded with 10 µL of internal standard solution. Calibration was performed using derivatized standard mixtures of varying concentrations added to the wells on the left side of the plate. Quality control samples were prepared by pooling equal volumes from each individual sample. All samples were stored at −80°C until analysis. Quantitative analysis of all target metabolites was conducted using an ultra-performance liquid chromatography–tandem mass spectrometry (UPLC-MS/MS) system (Acquity UPLC-Xevo TQ-S, Waters Corp., Milford, MA, USA). Raw data generated from the UPLC-MS/MS were processed using TMBQ software (v1.0, Metabo-Profile, Shanghai, China) for peak integration, calibration, and quantification of each metabolite. Metabolite concentrations in unknown samples were determined by comparison with a series of standards of known concentrations via calibration curves. These calibration curves characterize the relationship between the analytical signal and analyte concentration. For most metabolites, the instrument response (e.g., peak height or peak area) exhibits a linear relationship with concentration, which can be described by the equation y = ax + b, where y represents the instrument response, a is the slope reflecting the method's sensitivity, b accounts for the background signal, and x denotes the analyte concentration in the unknown sample. Using this model, metabolite concentrations in unknown samples were calculated from the measured instrument responses, enabling accurate and reproducible quantification.

In this study, we employed the Q300 kit (Metabo-Profile, Shanghai, China), an automated high-throughput metabolite array technology that enables quantitative detection of multiple metabolites across different concentration ranges within a single microtiter plate. This platform allows for absolute quantification as well as differential screening of diverse classes of metabolites, including amino acids, phenols, phenyl-derivatives or benzyl-derivatives, indoles, organic acids, fatty acids, carbohydrates, and bile acids. The kit incorporates isotope-labeled internal standards (e.g., L-arginine-15N2, hippuric acid-D5, TCDCA-D9, D-glucose-D7, carnitine-D3, C5:0-D9, and citric acid-D4), together with matched external standards, to ensure accurate qualitative and quantitative analyses (33). Of the more than 300 targeted metabolites, a total of 201 were successfully detected in our samples. The Q300 method provides broad coverage of disease-relevant metabolites, including those implicated in MDD and IBS, while offering high analytical accuracy, stability, and reproducibility.

All analyses in this study were performed within the R software environment. Clinical data were statistically analyzed using the R package tableone (v0.13.2) (34) . Categorical variables were analyzed using the chi-square test, and continuous variables were compared between the two groups using the t-test. Multiple hypothesis testing was adjusted using the Benjamini–Hochberg (BH) method to control the false discovery rate. Microbial community composition analysis was conducted using the R package vegan (v2.6-10) (35). Community structure differences were assessed by calculating Bray–Curtis distance matrices, and significance testing was performed via permutation-based multivariate analysis of variance (PERMANOVA) implemented by the adonis2 function. Alpha diversity indices were also calculated using the same package. Differential species analysis was conducted with the R package microeco (v1.9.1) in conjunction with the linear discriminant analysis effect size (LEfSe) method to identify characteristic taxa, with selection criteria of linear discriminant analysis (LDA) score >2 and P value <0.05 (36). Pathway enrichment analysis based on KEGG Orthology (KO) relative abundances was performed using the R package ReporterScore (v0.1.9), combining Wilcoxon and Kruskal–Wallis tests with 999 permutations. Significantly enriched pathways were selected based on absolute Reporter scores >1.96 and visualized accordingly (37). Correlations between variables were assessed using Spearman's rank correlation coefficients, with significance determined by corresponding P values.

To explore the co-variation among microbiome, serum metabolome, and psychological phenotypes, pairwise data sets were subjected to co-inertia analysis (CIA) using the R package made4 (v1.76.0) (38). After dimensionality reduction by principal component analysis (PCA), co-inertia structures were computed, and their significance was evaluated by permutation testing of the RV coefficient with 999 permutations. The results were visualized as arrow plots, illustrating sample trajectories across the two omics projection spaces, indicating synergistic variation. In the metabolomics analysis, the Wilcoxon rank-sum test was used to assess the significance between the two groups. PCA and orthogonal partial least squares discriminant analysis (OPLS-DA) were conducted for sample classification and modeling. Variable importance in projection (VIP) scores were subsequently calculated to evaluate the contribution of each metabolite to the model classification. Differential metabolites were identified based on the criteria of P value <0.05 and VIP > 1.

Among the 190 subjects, there were no statistically significant differences in age, sex, smoking status, alcohol consumption, BMI, religious affiliation, marital status, employment status, and monthly household income among the three groups (Table 1). The MDD group exhibited significantly higher scores on the IBS-SSS compared to the HC group. Furthermore, the MDD with IBS group demonstrated significantly elevated scores on the IBS-SSS as well as HAMD-17 and HAMA-14 scales when compared to the MDD group (Table 1). These findings suggest that our sampled MDD patients, those with comorbid IBS, had more pronounced levels of anxiety and depression.

One-hundred nineteen subjects were included for gut metagenomic sequencing and analysis, including 46 HC, 43 MDD, and 30 MDD with IBS. The total sequencing data for each sample and the proportion of reads classified as microbial are summarized in Table S1. Alpha diversity, assessed by Shannon and Simpson indices, was significantly higher in both MDD and MDD with IBS groups compared to HC. The Chao1 index was further elevated in the MDD with IBS group, indicating increased microbial richness under comorbid conditions (P < 0.05) (Fig. 2A). Principal coordinates analysis and PERMANOVA showed significant differences in microbial community structure between the disease groups and the HC group (P < 0.05), while no significant differences were observed between the MDD and MDD with IBS groups (Fig. 2B). Further PERMANOVA at the species level showed that group status significantly affected gut microbial beta diversity (R² = 0.0627, P < 0.001) while sex, age, BMI, smoking, and alcohol consumption showed no significant effects (P > 0.05) (Fig. 2C). At the phylum level, Firmicutes, Bacteroidota, and Actinobacteria were dominant across all samples, with Firmicutes most enriched in the disease groups (Fig. 2D). At the species level, the top 10 dominant species included Prevotella copri, Eubacterium rectale, and Phocaeicola vulgatus, with the HC group showing the highest total relative abundance of these species (Fig. 2E). LEfSe analysis identified numerous significantly different microbial taxa across taxonomic levels (from phylum to species) between groups (Fig. 2F). At the phylum level, Firmicutes were enriched in the MDD group, Bacteroidota in HC, and Actinobacteriota in the MDD with IBS group. At another level, several signature microbes were differentially enriched in each group. For instance, Prevotella copri was enriched in HC, while Clostridium scindens and Bifidobacterium animalis were enriched in the MDD with IBS group (Fig. 2G). Figure 2H illustrates the phylogenetic relationships of significantly enriched taxa with relatively high abundance across groups. Pairwise LEfSe analyses revealed significant microbial differences between disease groups and HC, but not between MDD and MDD with IBS (Table S2).

At the KO level, more differential KOs were identified between the MDD with IBS and HC groups (n = 63) than between the MDD and HC groups (n = 17). There were five shared upregulated KOs in both disease groups compared to HC (Fig. 3A). Most KOs uniquely altered in the MDD with IBS group also showed an upward trend in MDD, though to a lesser extent (Fig. 3B). KEGG pathway enrichment analysis revealed significant enrichment in multiple pathways in both MDD and MDD with IBS groups (|ReporterScore| > 1.96). Specifically, 47 pathways were enriched in MDD vs HC and 83 pathways in MDD with IBS vs HC (Fig. S1). In metabolism pathways, shared pathways primarily included those involved in carbon metabolism (map01200), pentose phosphate pathway (map00030), and 2-oxocarboxylic acid metabolism (map01210), among others. The MDD with IBS group uniquely enriched pathways, including D-amino acid metabolism (map00470), glycerolipid metabolism (map00561), and others (Fig. 3C and D). Further analysis showed stronger upregulation of functional genes in shared pathways in the MDD with IBS group. For example, only the multifunctional 2-oxoglutarate metabolism enzyme (K01616) was significantly upregulated in MDD, while eight KOs were upregulated in MDD with IBS (Fig. 3E and F). Notably, most enriched pathways were significantly correlated with emotional state and IBS-SSS score (Fig. S2). Although pathway enrichment analysis showed significant results between the two disease groups, no significant KO-level differences were observed, and no further pathway analysis was performed. These results suggest the necessity of integrating metabolomics to better understand the impact of gut microbiota on metabolic processes.

A total of 201 metabolites were detected across 119 serum samples (Table S3), and their classification is presented in Fig. 4A. PCA showed a significant difference between the MDD and HC groups (P = 0.044), whereas no significant differences were observed between the MDD with IBS and HC groups (P = 0.193) or between the MDD with IBS and MDD groups (P = 0.804) (Fig. 4B). Further analysis using an OPLS-DA model showed distinct clustering patterns among the groups (MDD with IBS vs HC: R²X = 0.136, R²Y = 0.62, Q² = 0.109; MDD with IBS vs MDD: R²X = 0.116, R²Y = 0.594, Q² = −0.647) (Fig. 4C). Based on the differential screening criteria (P < 0.05 and VIP > 1), multiple differential metabolites were identified. In the MDD with IBS group compared to the HC group, 10 metabolites were significantly upregulated and 19 were downregulated. Compared to the MDD group, only four metabolites were significantly downregulated in the MDD with IBS group (Fig. 4D). In addition, a total of nine metabolites showed significant differences exclusively in the MDD with IBS group.

We classified and analyzed the differential metabolites. The metabolite classification bar plot showed that most differential metabolites in the MDD with IBS group were derived from amino acids, fatty acids, carbohydrates, and bile acids (Fig. 4E). Analysis of the differential metabolites revealed that leucine, isoleucine, and valine were significantly downregulated in both the MDD vs HC and MDD with IBS vs HC comparisons. Meanwhile, glutaconic acid, glyceric acid, CDCA, GCDCA, and GCDCA-3S showed significant alterations only in the MDD with IBS vs HC comparison. Furthermore, four metabolites—N-acetylglutamine (NAG), imidazolepropionic acid, acetic acid, and N-acetylaspartic acid (NAA)—were significantly downregulated in the MDD with IBS vs MDD comparison. The box plots illustrate the relative abundance changes of these metabolites across the different groups (Fig. 4F through H).

Using CIA across all samples, we examined the influence of emotional state on gut microbiota composition and serum metabolite levels, as well as their interaction. Emotional state showed a marginal effect on both data sets, but the overall covariation between microbiota and metabolites was not statistically significant (Fig. 5A). MDD with IBS-specific differential microbes was identified, including 20 enriched and 10 depleted taxa, while the specific differential metabolites included only one enriched metabolite (glyceric acid) and seven depleted metabolites (Fig. 5B; Tables S2 and S4). Spearman correlation analysis revealed significant associations between several altered species and bile acid levels. Upregulated taxa such as Anaerotruncus colihominis, Lentihominibacter hominis, and Gordonibacter pamelaeae were positively correlated with CDCA and GCDCA levels. Additionally, Paraeggerthella hongkongensis and Blautia luti were negatively correlated with GCDCA. Conversely, downregulated species like Bacteroidetes bacterium ADurb Bin416 and Bacteroidales bacterium Barb4 showed positive correlations with GCDCA. No significant correlations were observed between glyceric acid and the altered microbiota (Fig. 5C). The abundances of bile acid–associated taxa are shown in Fig. 5D. Notably, shared differentially abundant taxa in both the MDD and MDD with IBS groups, such as Blautia obeum, Eggerthella lenta, and Clostridium scindens, are known to be involved in bile acid metabolism (Fig. 5E).

To further investigate the potential impact of microbial functional alterations on host metabolism, we performed integrative functional annotation and metabolite analysis in the MDD with IBS group. The differential metabolite glyceric acid was found to be involved in several key KO-enriched pathways, including the pentose phosphate pathway (map00030), glycerolipid metabolism (map00561), and carbon metabolism (map01200). The pentose phosphate pathway is a part of carbon metabolism, and glyceric acid is regulated by several enzymes, such as EC1.2.1.98, EC1.2.99.8, and EC1.2.7.5. The KO gene K03738, which encodes aldehyde ferredoxin oxidoreductase that catalyzes the EC1.2.7.5 reaction, was significantly upregulated in the MDD with IBS group (Fig. 6A). Taxonomic AOR (K03738) revealed that its primary microbial sources were the phyla Firmicutes and Actinobacteria, with some sequences only annotated at the other level, including Eggerthella and Oscillospiraceae. These taxa exhibited significant differential abundance in the MDD with IBS group (Fig. 6B). Further analysis of other key differentially expressed KO genes within this pathway (K08093 and K08094) indicated that they were predominantly derived from Firmicutes (Fig. 6C), notably from the families Eggerthellaceae, Enterococcaceae, and Lachnospiraceae, which also showed significant shifts in abundance (Fig. 6D). Collectively, these findings suggest that multiple functionally relevant genes implicated in glyceric acid-associated metabolic pathways originate from gut microbial taxa that differ significantly in abundance in MDD with IBS.

Lastly, an integrated network was constructed to explore the correlations among the top 20 significantly altered microbes, clinical indicators (emotional state and IBS-SSS scores), and key metabolites (including four bile acids and glyceric acid) in the MDD with IBS group. Several differential microbes showed significant correlations with bile acid levels and clinical scores, whereas glyceric acid exhibited generally weak correlations with these variables (Fig. 6E).