Participant characteristics

A total of 50 patients with depression and anxiety disorders were enrolled initially (Fig. 1b). After excluding participants who withdrew consent (n = 5), those who did not provide samples (n = 5), and those lost to follow-up (n = 6), 34 participants were included into the final analysis. Based on the HAM-A and HAM-D scores at time points 1 and 3, patients were designated as non-responders (n = 19), responders (n = 9), and remitters (n = 6). Nine of the 19 non-responders were further divided based on propensity score matching. The selected nine non-responders and nine responders were assigned to a discovery set, whereas the 10 non-responders who were excluded by propensity score matching, and the six remitters, were assigned to the validation set.

The clinical characteristics, including the dose of medication received during the observation period, of the participants are shown in Table 1. Factors used to calculate the propensity scores (age, sex, body mass index, and hospitalization (or not)) were matched successfully between propensity score-matched non-responders and responders. In the propensity score-matched discovery set, there was no significant difference between responders and non-responders with respect to treatment history during the study. The baseline HAM-A scores of non-responders were higher than those of responders. Since the validation set comprised non-responders and remitters, there were inevitable differences in the treatment strategies applied to each. The detailed clinical characteristics for each participant are shown in supplemental Table 1.

Dissimilarity among metabolomic and metagenomic data obtained from different participants was significantly higher than that among repeated samples obtained from the same participant, suggesting that the data were obtained appropriately (supplemental Fig. 1).

Table 1 Participant characteristicsElastic net analysis to discriminate between responders and non-responders in the discovery set

Elastic net analysis, used to compare baseline levels of metabolites between responders and non-responders in the propensity score-matched discovery set, identified 11 metabolites of interest (isethionate, N-ε-acetyllysine, cysteate, glycyl-L-leucine, taurine, guanosine, L-lysine, malonate, nicotinamide, indole-3-acetate, and dodecanedioate) (Fig. 2a, b, supplemental Table 2). In detail, predicted values based on the amount of each metabolite were calculated using the following formula: eY/(1 + eY)| Y = (-5.18 × 10− 3 × cysteate) + (-1.493 × 10− 4 × dodecanedioate) + (-4.063 × 10− 4 × glycyl-L-leucine) + (-1.29 × 10− 3 × guanosine) + (-7.083 × 10− 5 × indole-3-acetate) + (-7.11 × 10− 6 × isethionate) + (-4.487 × 10− 5 × L-lysine) + (-6.729 × 10− 4 × malonate) + (-3.065 × 10− 2 × N-ε-acetyllysine) + (2.569 × 10− 2 × nicotinamide) + (-4.665 × 10− 5 × taurine) + 1.308. The distribution of predictive values indicated that the combination of these metabolites could discriminate responders from non-responders in the propensity score-matched groups (Fig. 2c).

In the same way, elastic net analysis of metagenomic data identified 23 bacterial genera that differed in abundance between responders and non-responders in the propensity score-matched discovery set (Fig. 2d–e). Predicted values based on the amount of metabolites were calculated using the following formula: eY/(1 + eY)| Y = (26.7422 × f_Coriobacteriaceae_un) + (-3.0496 × Coprococcus) + (15220.1224 × k_Bacteria_Other) + (56.1757 × Corynebacterium) + (-0.2688 × Phascolarctobacterium) + (3854.053 × f_Peptostreptococcaceae_Other) + (1464.6944 × o_Streptophyta_un) + (13.8465 × Parabacteroides) + (660.2579 × Oxalobacter) + (1.7989 × Prevotella) + (-120.8179 × cc_115) + (-138.6197 × rc4-4) + (-735.3305 × SMB53) + (3409.7147 × Atopobium) + (14.3232 × f_S24-7;g_unknown) + (7896.5654 × f_Streptococcaceae_Other) + (-11.5373 × f_Clostridiaceae_un) + (14392.5596 × Abiotrophia) + (-2196.5849 × o_Lactobacillales_un) + (4.2996 × Desulfovibrio) + (-1115.5825 × f_Clostridiaceae_Other) + (3.3706 × Odoribacter)+ (-5.0804 × Clostridium) -0.106. The distribution of predictive values indicated that the combination of these genera could discriminate responders from non-responders in the propensity score-matched groups (Fig. 2f). After removing these genera, elastic net analysis identified a single genus, Bilophila, that could discriminate responders from non-responders (AUROC = 0.88).

Fig. 2

Variable selection using elastic net analysis. (a) The coefficient path for the elastic net-regularized logistic regression applied to the metabolome data dependent on log(?). The number of non-zero coefficients is shown above the plot. (b) Cross-validation binomial deviance curve derived from elastic net-regularized logistic regression analysis of metabolome data, along with one-standard-error bands calculated from 10-fold realizations. The vertical line corresponds to the minimum value for log(?). (c) Violin plot of the predicted values calculated by elastic net-regularized logistic regression analysis of the metabolome data. (d) Coefficient path for the elastic net-regularized logistic regression applied to the metagenome data dependent on log(?). The number of non-zero coefficients is shown above the plot. (e) Cross-validation binomial deviance curve for the elastic net regularized logistic regression analysis of metagenome data, with one-standard-error bands calculated from 10-fold realizations. The vertical line corresponds to the minimum value for log(?). (f) Violin plot of the predicted values calculated by elastic net-regularized logistic regression analysis of the metagenome data

Identification of metabolites associated with response to treatment

ROC analysis of the 11 identified metabolites revealed that the baseline levels of nine (isethionate, N-ε-acetyllysine, cysteate, glycyl-L-leucine, taurine, guanosine, L-lysine, malonate, and nicotinamide) could be used to distinguish between responders and non-responders with statistical significance (Fig. 3a). A comparison of the average baseline levels allowed us to narrow this down further to seven metabolites (N-ε-acetyllysine, cysteate, glycyl-L-leucine, taurine, guanosine, L-lysine, and malonate) (Fig. 3b), the levels of which were lower in responders than in non-responders.

Fig. 3

Identification of metabolites associated with therapeutic efficacy. (a) Receiver operating characteristic (ROC) analysis of the identified fecal metabolites. The area under the curve (AUC) and 95% confidence intervals (CIs) are shown for each curve. Data shown in red are statistically significant. Yellow-green areas denote the 95% CIs for the ROC curves. (b) Violin plots of the baseline levels of identified fecal metabolites. P-values for responders versus non-responders were calculated using Welch’s t-test. Data shown in red are statistically significant

Correlation between gut microbes and the identified metabolites

To infer the origin of the identified metabolites, we used the genera identified by elastic net analysis of the propensity score-matched discovery set to compare the baseline microbial composition in the gut of responders and non-responders (Fig. 4a). At the phylum level, Bacillota were dominant in non-responders, while Bacteroidota were dominant in responders. This suggests that the low levels of the seven identified metabolites may be associated with a reduction in Bacillota.

We then focused on 13 genera (Coprococcus, Prevotella, Parabacteroides, f_Clostridiaceae_un, Phascolarctobacterium, Odoribacter, f_Coriobacteriaceae_un, Clostridium, f_S24-7_un, Desulfovibrio, cc_115, Corynebacterium, and Bilophila) that accounted for a relatively large proportion among every identified genus (supplemental Fig. 2). Then, we investigated correlations between the identified metabolites and genera in the discovery set at time point 1 (Fig. 4b). The levels of L-lysine and N-ε-acetyllysine correlated positively with the abundance of Coprococcus (R = 0.56 [95% CI, 0.13–0.82] and 0.52 [95% CI, 0.07–0.79], respectively); the levels of cysteate correlated positively with the abundance of f_Clostridiaceae_un (R = 0.51 [95% CI, 0.19–0.84]); and the levels of taurine correlated positively with the abundance of Phascolarctobacterium (R = 0.53 [95% CI, 0.09–0.80]). All three of these genera belong to Bacillota. By contrast, the levels of L-lysine, glycyl-L-leucine, and N-ε-acetyllysine were correlated negatively with the abundance of Odoribacter (R = -0.55 [95% CI, -0.81– -0.11], -0.64 [95% CI, -0.85– -0.25], and − 0.54 [95% CI, -0.81– -0.10], respectively), which belong to Bacteroidota. The levels of cysteate and taurine correlated negatively with the abundance of Bilophila (R = -0.66 [95% CI, -0.86– -0.27] and − 0.59 [95% CI, -0.83– -0.18], respectively), which belong to Pseudomonadota, whereas the levels of L-lysine were correlated negatively with the abundance of f_Coriobacteriaceae_un (R = -0.50 [95% CI, -0.78– -0.04]), which belong to Actinomycetota.

Subsequently, we checked whether the same correlations between genera and metabolites were true for the validation set (Fig. 4c). Although the correlations between all investigated pairs did not reach statistical significance, the correlation coefficient for the levels of Odoribacter and N-ε-acetyllysine (R = -0.47 [95% CI, -0.78–0.03]) was similar to that of the discovery set.

Since the baseline HAM-A scores between responders and non-responders were significantly different, we compared the ability of baseline HAM-A scores, N-ε-acetyllysine levels, and the abundance of Odoribacter to predict responses to treatment. The AUROC for baseline N-ε-acetyllysine levels (Fig. 3a) was greater than that for baseline HAM-A scores or that for the abundance of Odoribacter (supplemental Fig. 3), suggesting that measuring N-ε-acetyllysine levels can assist prediction of treatment responses.

Fig. 4

Association between gut microbes and identified metabolites. (a) Differences in the composition of the gut microbiota between responders and non-responders. Genera with high bacterial abundance are shown in red (See supplemental Fig. 2). (b) Correlations between levels of the identified metabolites and the abundance of the identified genera. The numbers in the boxes are Pearson’s correlation coefficients. Yellow highlights denote statistical significance. (c) Correlation between metabolite levels and the abundance of genera in the validation set. Pearson’s correlation coefficients (R) and 95% confidence intervals (CIs) are shown for each curve. Data shown in red indicate |R|>0.4. Yellow-green areas denote the 95% CIs for the regression lines. (d) Changes in N-ε-acetyllysine levels, and in the abundance of Odoribacter, over time. P-values, calculated using the linear mixed model with Bonferroni correction, denote the significance of differences in trends between groups

Changes in metabolite levels and in the abundance of Odoribacter during the study

We examined changes in the levels of the identified metabolites in responders, non-responders, and remitters during the study period. The levels of N-ε-acetyllysine were significantly lower in responders than in the other groups at time point 1, whereas the levels in responders increased to the same levels as those in non-responders and remitters at time point 3 (Fig. 4d), suggesting that recovery of N-ε-acetyllysine levels may be related to an improvement in symptoms. The levels of L-lysine showed a similar trend (supplemental Fig. 4). By contrast, the abundance of Odoribacter was higher in responders than the other groups at time point 1. This difference was sustained throughout the study period, suggesting that there is no clear association between a change of bacterial composition and improvement of symptoms in patients with depression and anxiety.

Finally, we examined changes in the levels of the identified metabolites, and in the abundance of Odoribacter, under different treatment strategies (e.g., the use of antipsychotics, hospitalization, and the electroconvulsive therapy [ECT]) during the study period (Fig. 5). Interestingly, N-ε-acetyllysine levels were elevated significantly in participants who were hospitalized and underwent ECT. Antipsychotic use did not affect N-ε-acetyllysine levels. Antidepressants were used by almost all participants, so a meaningful analysis of their effects was not possible. Thus, changes in intestinal metabolites may be caused by changes in lifestyle, such as dietary content due to hospitalization. Measurement of N-ε-acetyllysine levels may also help to determine whether inpatient treatment is required.

Fig. 5

Effects of treatment strategy on metabolites and gut microbes. Changes in N-ε-acetyllysine and L-lysine levels, and in the abundance of Odoribacter, over time in patients receiving different treatments. P-values, calculated using the linear mixed model, indicate the significance of differences in trends between groups