Antidepressants affect gut microbiota composition

To investigate whether antidepressants may alter gut microbiota, we chose five different antidepressants common in clinical practice and different in their mode of action. We used two selective serotonin reuptake inhibitors (SSRIs) - fluoxetine and escitalopram, two serotonin norepinephrine reuptake inhibitors (SNRIs) - venlafaxine and duloxetine, and desipramine that acts as norepinephrine reuptake inhibitor. Since antidepressants require at least three weeks to show their therapeutic effects, we characterized the gut microbial community after 21 day of daily antidepressant treatment (Fig. 1a; OTU table with all detected bacteria can be found in Supplementary Table 1). Antidepressants were injected i.p. in order to provide a specific concentration of antidepressants directly to the gut. Treatment by drinking water would cause potential cofounders in the analysis of specific drug effects on microbiota, due to potential variability among experimental groups in drinking water volumes and differential kinetics of breakdown of antidepressants in drinking water. Nonetheless, one important limitation of our approach is that daily i.p. administration adds an element of stress to the experimental outline which can affect experimental findings.

Analyses of alpha diversity revealed that all antidepressants, except desipramine, reduced the richness of microbial communities (Fig. 1b; Supplemental Figure S1a, b), but did not affect their evenness (Fig. 1c). Beta diversity measures of gut microbiota were also affected by all studied antidepressants, but more pronounced effects were observed in unweighted UniFrac analyses (Fig. 1d–f) than in weighted UniFrac (Supplemental Figure S1c–e). Namely, beta diversity of fecal microbial communities from mice receiving antidepressants was higher than beta diversity of control samples (Fig. 1d, f). Further, microbial communities of control group were more similar to each other than when they were compared to samples of any of antidepressant treated groups (Fig. 1e, f).

In order to identify bacterial taxa which differed in relative abundances in antidepressant treated mice in comparison to controls, we analyzed the 16S rRNA gene sequencing results using linear discriminant analysis (LDA) effect size (LEfSe) algorithm. A comparison between control mice and all antidepressant groups together revealed that Ruminococcus, Adlercreutzia and an undefined genus in the order RF32, class Alphaproteobacteria were less abundant in antidepressant treated mice (Fig. 2a, b). When analyzing the effects of each individual antidepressant in comparison to control using pairwise comparisons, the same genera were shown to be less abundant in escitalopram, venlafaxine, duloxetine and desipramine groups, but not in fluoxetine group (Supplemental Figure S2). Further, in both the overall effects of antidepressants and in the pairwise analysis, the differences in genus Adlercreutzia contributed to observed decreased abundance of family Coriobacteriaceae, order Coriobacteriales, class Coriobacteriia, phylum Actinobacteria in antidepressant groups (Fig. 2a, b; Supplemental Figure S2). In addition to the LEfSe analysis, we further analyzed the 16S rRNA gene sequencing data with PERMANOVA (Supplementary Table 2). At the genus level, both Ruminococcus and Adlercreutzia were effected by antidepressant treatment (p < 0.05), although only Adlercreutzia remained significant after corrections for multiple comparisons. In direct pairwise analysis between control and each of the antidepressant groups, some antidepressants significantly decreased levels of Ruminococcus, although they were no longer significant after corrections for multiple comparisons. Therefore, the findings of decrease in Ruminococcus and Adlercreutzia were more significant in the LEfSe analysis, and we further explored these findings using qRT-PCR.

Validation of antidepressant-induced decrease in Ruminococcus and Adlercreutzia levels was done by qRT-PCR (undefined genus in order RF32 could not be analyzed by qRT-PCR because of the shortage of knowledge about its sequence). In order to determine which species of Ruminococcus was most affected, examination of the OTU table revealed that the abundance of the OTU assigned to Ruminococcus flavefaciens was altered by most of the antidepressant treatments (Kruskal–Wallis test H = 27.88, p < 0.001) (Fig. 2c). The qRT-PCR verified that R. flavefaciens was indeed less abundant in all antidepressant groups compared to controls (Kruskal–Wallis test H = 22.33, p < 0.001) (Fig. 2d). In the genus Adlercreutzia, Adlercreutzia equolifaciens is the only characterized species of this genus. The OTU that had the highest (98%) similarity to Adlercreutzia equolifaciens and was the most abundant in mice gut, was also decreased by most of antidepressants (Kruskal–Wallis test H = 16.06, p < 0.01) (Fig. 2e). Likewise, the qRT-PCR results showed that all antidepressants except desipramine, reduced levels of A. equolifaciens in the mice gut (Kruskal–Wallis test H = 26.60, p < 0.001) (Fig. 2f).

R. flavefaciens but not A. equolifaciens reduces antidepressive effects of duloxetine

In our next group of experiments, we wanted to explore the hypothesis that R. flavefaciens or A. equolifacien may mediate the effects of antidepressant treatment on depressive-like behavior. With that aim, we treated BALB/c mice with antidepressant, bacteria (R. flavefaciens or A. equolifaciens), or both, and then performed behavioral testing (Fig. 3a). For this group of experiments, we chose antidepressant duloxetine because it decreased both R. flavefaciens and A. equolifaciens, as well as it showed the biggest effect in TST in BALB/c mice (Supplemental Figure S3). Presence of gavaged bacteria in the gut was confirmed by detecting R. flavefaciens or A. equolifaciens in stool samples of treated mice (Supplemental Figure S4).

First, we examined effects of R. flavefaciens in the TST. As expected, duloxetine induced a significant effect (two-way ANOVA: F dul = 30.70, p < 0.001), and duloxetine treated animals displayed significantly less immobility, compared to control animals (Fig. 3b). Interestingly, there was also a significant effect of R. flavefaciens (two-way ANOVA: F Rum = 5.90, p < 0.05), and animals treated with both R. flavefaciens and duloxetine displayed significantly more immobility compared to those treated with duloxetine alone (Fig. 3b). Therefore, R. flavefaciens treatment was able attenuate the duloxetine effect in the TST. In the FST, significant effects were also obtained by both duloxetine (two-way ANOVA: F dul = 5.70, p < 0.05) and R. flavefaciens (two-way ANOVA: F Rum = 8.74, p < 0.01) (Fig. 3c). Specifically, animals concomitantly treated with R. flavefaciens and duloxetine were more immobile than animals treated with duloxetine alone, and displayed behavior comparable to control animals (Fig. 3c). Therefore, R. flavefaciens treatment was able to abolish the effect of duloxetine treatment in the FST. The effects of R. flavefaciens in TST and FST were not confounded by locomotor deficits (Supplemental Figure S5a, b). Together, these results demonstrate that R. flavefaciens can attenuate duloxetine effects in behavior despair paradigms.

Further, we examined the effect of R. flavefaciens on anhedonia, by SPT. Both the antidepressant and the bacteria showed significant effects (two-way ANOVA: F dul = 6.42, p < 0.05; F Rum = 13.73, p = 0.001) (Fig. 3d). Namely, duloxetine increased mice preference for 2% sucrose, while it did not change the sucrose preference in the group receiving R. flavefaciens together with the drug (Fig. 3d). Also, sucrose preference was higher in group receiving duloxetine than in group receiving both duloxetine and R. flavefaciens (Fig. 3d). Moreover, the interaction of duloxetine and the bacteria treatment was significant in SPT (two-way ANOVA: F dul*Rum = 8.34, p < 0.01) (Fig. 3d). In conclusion, R. flavefaciens supplementation reduced or abolished antidepressive properties of duloxetine.

In addition, the effect of R. flavefaciens on general gastrointestinal health was evaluated by number of fecal pellets in open field test (Supplemental Figure S5c). R. flavefaciens abolished constipation induced by duloxetine. The decreased defecation in a new environment induced by duloxetine can be attributed to its anxiolytic effect as well, yet we note that there were no changes in time spent in center of the open field arena by any of treatments, which is also taken as sign of anxiety behavior54 (data not shown).

Additionally, we determined whether antidepressant effects can be modulated by A. equolifaciens. As in the previous experiment, duloxetine had significant effects on immobility time in TST (two-way ANOVA: F dul = 55.63, p < 0.001) and to lesser degree in FST (two-way ANOVA: F dul = 6.52, p < 0.05) (Fig. 3e, f). However, when the drug was given together with A. equolifaciens, it still exhibited its antidepressant properties (Fig. 3e, f). Therefore, contrary to R. flavefaciens, A. equolifaciens did not abolish antidepressive effects of duloxetine.

R. flavefaciens up-regulates mitochondrial genes while down-regulates neural genes in mPFC

In order to reveal mechanisms by which R. flavefaciens exerts its effects on the brain, we did whole transcriptome analyses on RNA extracted from the mPFC from each of the four experimental groups: control, duloxetine treated, R. flavefaciens treated, and those treated with both R. flavefaciens and duloxetine. The mPFC is known to be a significant center of behavior regulation, receiving inputs from different limbic structures, and has been implicated in both depression and antidepressant treatment effects55,56.

Differential expression analyses

RNA-seq data were firstly analyzed by differential expression analyses (DEA) (Supplementary Table 3). Duloxetine treatment alone, compared to control, changed expression of only one gene, Adrb1, one of the norepinephrine receptors (Fig. 4a). In contrast, R. flavefaciens treatment resulted in 324 differentially expressed genes (DEGs), in comparison to controls (Fig. 4a). GO analyses revealed that genes up-regulated by R. flavefaciens treatment were enriched for mitochondrial processes (Fig. 4b), especially oxidative phosphorylation, while the down-regulated genes were enriched with genes involved in neural plasticity (Fig. 4c). When the group treated by both R. flavefaciens and duloxetine was compared to controls, there were 185 DEGs (Fig. 4a). Interestingly, concomitant R. flavefaciens and duloxetine treatments also resulted in the up-regulation of genes enriched for mitochondrial energy metabolism (Fig. 4d), and down-regulation of genes enriched for neural plasticity (Fig. 4e). When the group receiving both the bacteria and duloxetine was compared to the group receiving duloxetine alone, only one gene, Dpysl2, was shown to be differentially expressed, i.e., down-regulated. The list of all DEGs can be found in Supplemental Table 3. Therefore, treatment with R. flavefaciens induces expression changes of genes related to mitochondrial and neuronal process in the mPFC.

Fig. 4: Differentially expressed genes (DEGs) after duloxetine and R. flavefaciens treatments. a The total number of DEGs, as well as the number of genes that are up-regulated and down-regulated by the particular treatment, are represented in the table. b, c Gene ontology (GO) enrichment analyses of genes up-regulated (b) and down-regulated (c) by R. flavefaciens treatment. d, e GO enrichment analyses of genes up-regulated (d) and down-regulated (E) by concomitant duloxetine and R. flavefaciens treatment. Bars representing GO terms show Benjamini and Hochberg FDR adjusted p values. n = 6 animals per experimental group. dul duloxetine, Rum R. flavefaciens Full size image

Weighted gene correlation network analysis

In order to determine clusters of genes whose expression was associated with duloxetine and R. flavefaciens treatment, we analyzed our data using weighted gene correlation network analysis (WGCNA). Namely, DEA are useful method for independent discovery of genes, with high confidence, which expression levels are changed by a treatment. On the other hand, WGCNA approach enabled us to discover groups of genes, which individually may not have been identified by DEA, but whose changes in expression are strongly mutually correlated and related with the same direction to a treatment. Further, this method could be of greater relevance for understanding disturbed biological pathways.

WGCNA clustering of our RNA-seq data identified a total of 18 modules of co-expressed genes. Module-trait relationship revealed that there were three modules significantly correlated with R. flavefaciens treatment that were not correlated with the same directionality with duloxetine treatment (Fig. 5a). The blue module (2904 genes), was positively correlated to the bacterial treatment, while the turquoise module (3353 genes), as well as a small skyblue module (186 genes), were negatively correlated to the bacterial treatment (Fig. 5a). Further confirmations of R. flavefaciens effects were strong correlations of module membership (MM) with gene-trait relationship for all three modules (Supplemental Figure S6). In other words, the genes that were strongly correlated to the module eigengenes (i.e., genes with high MM) were also strongly correlated with presence/absence of R. flavefaciens treatment.

Fig. 5: Results of weighted gene correlation network analysis (WGCNA). a Table of module-trait relationship reports Kendall’s correlation coefficients, and its corresponding p values, between the eigengene value of each module and the particular treatments. Modules related to R. flavefaciens treatment, and were not correlated with the same directionality with duloxetine treatment, are emphasized by dashed lines. b–g Gene ontology (GO) and protein-protein interaction (PPI) network analysis of WGCNA modules significantly related to R. flavefaciens treatment. GO enrichment analyses (b) and PPI network analysis (c, d) of genes in blue module with module membership (MM) > 0.7. GO enrichment analyses (e) and PPI network analysis (f) of genes in turquoise module with MM > 0.7. GO enrichment analyses of genes in skyblue module, with MM > 0.7 (g) Bars representing GO terms show Benjamini and Hochberg FDR adjusted p values. Node size in PPI networks indicates number of interactions with other nodes in the network (i.e., degree centrality), while the node color reflects number of shortest paths that rely on that given node within the network (i.e., betweenness centrality). dul duloxetine, Rum R. flavefaciens Full size image

In order to understand biological mechanisms associated with genes in these three modules, we performed GO and protein-protein interaction (PPI) network analyses on module genes (MM > 0.7). GO analyses of genes in the blue module, which was positively correlated with R. flavefaciens treatment, determined enrichment of mitochondrial energy generation processes (Fig. 5b). Subsequent PPI network analyses of the blue module revealed two protein clusters of strongly interconnected gene products also involved in mitochondrial functions. One protein cluster (composed of 32 nodes) mainly consisted of subunits belonging to mitochondrial complex I (NADH:ubiquinone oxidoreductase (NDUF)), while 4 genes are belonging to mitochondrial complex 3 (ubiquinol-cytochrome c reductase complex subunits (UQCR)) (Fig. 5c). The members of the other protein cluster (composed of 14 nodes) are part of mitochondrial ribosomal protein (MRPL) gene family, involved in translation of mitochondrial genes, that are subunits of respiratory chain complex (Fig. 5d). Therefore, R. flavefaciens treatment was associated with the upregulation of several mitochondrial pathways.

Turquoise and skyblue modules, whose expression levels were negatively related to R. flavefaciens treatment, were both enriched with genes involved in neural functions (Fig. 5e, g). PPI analysis of turquoise module revealed one cluster of 51 proteins, consisting of several smaller highly interacting clusters (Fig. 5f). Gene products of these protein clusters are involved in ionotropic glutamate neurotransmission (such as Gria1 and Gria2, AMPA receptor subunits, and Grin2a and Grin2b, NMDA receptor subunits), protein phosphorylation (such as Prkacb, a subunit of protein kinase A (PKA)), phosphatidylinositol signaling (such as Pten) and vesicle-mediated trafficking (such as Dync1h1 and Uso1). Genes from the smallest skyblue module did not give any relevant protein network. Therefore, the WGCNA provided further evidence that R. flavefaciens treatment affects gene expression pathways involved in neuronal function, and further highlight the glutamatergic system.

Modules affected by duloxetine treatment (Fig. 5a) were enriched for genes involved in RNA splicing (paleturquise module), transcriptional regulation and MAPK activity (lightgreen module) and protein ubiquitination and cellular localization (tan module) (Supplemental Figure S7), while cyan and darkgrey module did not show enrichment for any GO terms. Therefore, duloxetine also had interesting effects on gene networks. However, these networks were not affected by R. flavefaciens, and do not appear to be the target of R. flavefaciens effects.

Together, bioinformatical analyses of mPFC transcriptome showed that R. flavefaciens treatment up-regulated genes involved in mitochondrial oxidative phosphorylation, while down-regulating genes involved in neural plasticity. Therefore, these putative pathways may explain the modulatory effects of R. flavefaciens on antidepressant actions on depressive-like behavior.

Effects of R. flavefaciens on serotonin and noradrenalin effects in mPFC