Cholic acid (CA), chenodeoxycholic acid (CDCA), ursodeoxycholic acid (UDCA), deoxycholic acid (DCA), lithocholic acid (LCA), taurocholic acid (TCA), taurochenodeoxycholic acid (TCDCA), tauroursodeoxycholic acid (TUDCA), taurodeoxycholic acid (TDCA), taurolithocholic acid (TLCA), glycocholic acid (GCA), glycochenodeoxycholic acid (GCDCA), glycoursodeoxycholic acid (GUDCA), glycodeoxycholic acid (GDCA), and lithocholic acid-2,2,3,4,4-d5 (internal standard) were all purchased from Sigma-Aldrich (St. Louis, MO, USA). Glycolithocholic acid (GLCA) was purchased from J&K Scientific Ltd. (Shanghai, China). HPLC-grade ammonium formate (≥99%), ammonium acetate, methanol, and acetonitrile were purchased from Merck KGaA (Darmstadt, Germany). HPLC-grade formic acid (99%) was purchased from Anaqua Chemicals Supply (Wilmington, USA).

Seven-week-old male ICR mice were purchased from SLAC Laboratory Animal Co., Ltd. (Shanghai, China). After their arrival, mice were single-caged and divided into the normal control group and the CUMS model group of 12 animals each randomly based on their body weight and sucrose preference test results. Mice were acclimated for 7 days in a temperature- (23–26°C) and humidity-controlled (40–60%) room under a 12-h light/dark cycle (lights on 07:00–19:00) with free access to food and water. During the modeling period, mice were weighed biweekly.

CUMS progression contained a total of 8 different stimulations including: 1) food deprivation for 24 h, 2) water deprivation for 24 h, 3) damp sawdust for 24 h, 4) tail pinching for 2 min, 5) restraint for 1 h, 6) cage tilting at 45° for 24 h, 7) cold swimming for 10 min, and 8) day and night reversal for 24 h. Two or three types of stimulations were delivered daily and randomly to the mice in the model groups for 56 days.

Depression-related behavioral tests including the sucrose preference test (SPT), forced swim test (FST), and tail suspension test (TST) were performed during the experimental period.

For the SPT, all mice were habituated to 1% sucrose solution during the adaptation cycles. After the adaptation progression, mice were deprived of water and food for 12 h and were provided with free access to two tubes containing 20 ml of sucrose solution (1% w/v) and water respectively for 5 h. The sucrose preference rate was calculated subsequently using the following formula: sucrose preference = volume of sucrose consumed/total volume (water and sucrose) consumed × 100%. We performed the SPT on day 57 to evaluate the modeling effect.

We conducted the forced swim and tail suspension tests on days 58 and 59, respectively. During the forced swim test, the mice were individually placed into glass cylinders (height of 40 cm, diameter of 18 cm) containing 25°C water at a depth of 15 cm for 10 min. Immobility time was measured of last 4 min was recorded to estimate the symptom of depression. During the tail suspension test, mice were individually suspended by their tails for 6 min using a small piece of tape on the shelf, placed at the height of 60 cm above the floor. The duration of immobility during the final 4 min was recorded to measure depressive status.

On day 60, serum and feces were collected after 12 h of fasting. Livers and intestinal contents were removed immediately after the mice had been sacrificed. All samples were stored in a freezer at −80°C for further processing.

Serum samples were thawed on ice and 400 µl of methanol was subsequently added into 100 µl of serum sample in an EP tube. The mixture was vortexed for 1 min and centrifuged at 15,000 g for 10 min at 4°C. The supernatant was then transferred into another EP tube and evaporated to dry with an Eppendorf Vacufuge Concentrator 5305. The residue was resuspended in 150 ul of 80% methanol and filtered through a 0.22-µm nylon syringe filter. For all samples, equal volumes of solutions were mixed into quality controlled samples to evaluate instrument analysis stability and repeatability.

The separation of the target compounds was performed on a Waters ACQUITY UPLC HSS T3 (2.1 mm × 150 mm, 1.8 µm) liquid chromatography column at 40°C with a ACQUITY UPLC CSH C18 VanGuard Pre-column (2.1 mm × 5 mm, 1.7 µm) using a Dionex Ultimate 3000 UPLC system. The mobile phase contained 0.1% aqueous formic acid and 0.1% formic acid in acetonitrile in positive ion mode and 5 mM ammonium formate aqueous buffer and acetonitrile in negative ion mode. The mobile phase flow rate was 0.25 ml/min and the injection volume was 5 µL both in the positive and negative ion modes. Supplementary Table S1 presents the detailed gradient elution conditions. The Q Exactive Orbitrap mass spectrometer (Thermo Fisher Scientific, USA) equipped with an ESI interface was applied for mass spectrometry analysis. The optimal parameters were as follows: sheath gas flow rate, 30 arb; aux gas flow rate, 10 arb; capillary temperature, 325°C; scan range: 81–1000 Da; stepped normalized collision energy, 30 in NCE mode; spray voltage, 3.5 kV (positive)/−2.5 kV (negative). All MS spectra were acquired and analyzed using the Xcalibur 4.0 software (Thermo Fisher Scientific).

Metabolomics analysis was carried out by BioNovoGene (Suzhou, China). After the raw data files were converted into an mzXML format by the ProteoWizard software (v3.0.8789), the freely available XCMS software was used to perform peak identification, filtration, alignment, and integration. The three-dimensional data matrix, including retention time, mass to charge ratio, and intensity, was converted into a table for further process analysis. In order to compare the data of different magnitudes, the peak area of the data was batch-normalized before multivariate statistical analysis. The data were then uploaded into SIMCA-P 13.0 to perform principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA). Autoscaling was used in all the models to achieve more scientific, reliable, and intuitive results. The variable importance in the project values (VIP) obtained from the OPLS-DA model and p-value from Student’s t-test were used to select the potential metabolites in the study. Metabolites with VIP >1 and p-value < 0.05 were considered statistically significant. These potential metabolites were subsequently subjected to pathway analysis performed through MetaboAnalyst 5.0 (https://www.metaboanalyst.ca/).

The bile acid profiles in the serum, feces, and liver were quantified using our previously validated UPLC-Q/Orbitrap-HRMS methods. Briefly, the bile acids in the feces and liver were extracted with 5 vol of deionized water by Qiagen TissueLyserII. For the feces samples, 200 uL of acetonitrile and 10 uL of 14% ammonia solution were added into 100 uL of fecal suspensions spiked with internal standard. For the serum and liver homogenates, 400 uL of acetonitrile was added into 100 uL of sample spiked with internal standard. The mixture was vortexed for 1 min and centrifuged at 20,000 g for 10 min at 4°C. The supernatant was then transferred into another EP tube and evaporated to dry with a vacuum centrifugal concentrator. The residue was resuspended in 100 uL of 80% methanol and filtered through a 0.22-µm nylon syringe filter. All blank matrix used for the calibration standard configurations and quality control samples were prepared using the activated carbon adsorption method.

The separation of the target compounds was performed on the same instrument and column as described above. The mobile phase flow rate and injection volume were 0.2 ml/min and 5 μL, respectively, with 10 mM ammonium formate aqueous buffer (A) and acetonitrile (B). The optimized gradient elution (0–7 min, 35–60% B; 7–8.5 min, 60–95% B; 8.5–12 min, 95% B; 12–12.3 min, 95%–35% B; 12.3–16 min, 35% B) was performed to separate the different bile acid components. Acquisition was performed in negative selective ion monitoring mode. All MS spectra were acquired and analyzed using the Xcalibur 4.0 software.

Total genomic DNA from the intestinal contents (5 samples from each of the control and model group) was extracted using the HiPure Stool DNA Kit (Megan, Guangzhou, China) according to the manufacturer’s protocols. The DNA concentration was measured using the Equalbit dsDNA HS Assay Kit (Novizan, Nanjing China). The NGS library preparation and Illumina sequencing was performed by GENEWIZ, Inc. (Suzhou, China). Approximately 20–30 ng of DNA was used to generate amplicons. The V3 and V4 hypervariable microbial 16S rDNA regions were amplified by PCR using a panel of proprietary primers designed by GENEWIZ. A linker with an index was then added to the end of the PCR product of 16S rDNA by PCR for NGS sequencing. The obtained sequencing library was subsequently purified with magnetic beads, followed by library quality control checks using a microplate reader and agarose gel electrophoresis. The library was then quantified to 10 nM and PE250/FE300 paired-end sequencing was performed using an Illumina MiSeq instrument (Illumina, San Diego, CA, USA).

Next, the forward and reverse reads were joined in pairs, followed by filtering the sequences containing N in the splicing results and retaining the sequences with a length beyond 200 bp. The obtained longer sequences were used to perform sequence clustering using VSEARCH (1.9.6) (sequence similarity was set to 97%) against the Silva_138 16SrRNA database (http://www.arb-silva.de/). The Ribosomal Database Program classifier was used to assign taxonomic categories to predict the community composition at the genus levels. Sequence data associated with this project have been deposited in the NCBI database (Accession Number: PRJNA796629).

All statistical analyses were performed using the GraphPad Prism 9.0 software. A two-tailed t-test was performed to compare between the groups and statistically significant differences were labeled with one, two, three, or four asterisks corresponding to p < 0.05, p < 0.01, or p < 0.001, respectively. Correlations between the gut microbiotic abundance and bile acid profiles were estimated using Pearson’s correlation analysis.

The body weight of the animals was measured before and during the treatment period. The mice in the model group gained less weight than control group at the end of the CUMS progression (Figures 1A–D). In addition, significant CUMS effects were present in the case of the sucrose consumption in SPT, immobility time in both FST and TST compared with the control group. The results demonstrated that a CUMS mouse model was successfully created.

To investigate the impact of chronic stress upon the metabolomic profiling, untargeted metabolomic analysis was performed to analyze the metabolite composition in the serum of mice. The obvious separation trend from the PCA (Supplementary Figure S1) and the OPLS-DA (Figures 2A,B) score plot indicated metabolic differences between the groups. Our OPLS-DA permutation test showed that the model we established was not over-fitting (Figures 2C,D).

In order to screen out potential metabolites, we used the VIP value of the OPLS-DA model beyond 1.0 and the p-value of the two-tailed unpaired Student’s t-test results less than 0.05 as a threshold to distinguish the metabolites from the model and control groups. A total of 74 metabolites were significantly changed during the CUMS progression. Table 1 shows the detailed information of these metabolites. The affected pathways mainly involved amino acid, sugar, nucleotide metabolism, unsaturated fatty acid biosynthesis and metabolism, vitamin synthesis and absorption, and bile acid metabolism. Supplementary Figure S2A and Supplementary Table S2 show the bubble chart of the KEGG pathway analysis and the detailed information of the pathway analysis, respectively. We found that the two main primary bile acid (chenodeoxycholic acid and taurocholic acid) levels were significantly altered in the model group, indicating abnormalities in bile acid synthesis or metabolism.

To further examine the bile acid metabolism disrupted by CUMS progression, we quantified the detailed bile acid profiles in the serum using our previously established method (Supplementary Figure S3 shows the chromatographic separation of the different components). CUMS significantly increased the level of three free bile acids, UDCA (345%↑, p = 0.0314), CDCA (220%↑, p = 0.0152), and DCA (197%↑, p = 0.0009), whereas it significantly reduced the level of taurine-conjugated primary bile acid TCA (56%↓, p = 0.0452) (Figure 3A). The taurine-conjugated-to-free bile acid ratios in the model group were significantly lower than those in the control group (Figure 4A). In addition, the hydrophobicity index (HI) of the circulating bile acid pool was calculated as described previously (Heuman, 1989). We observed that the HI in the model group was significantly higher than that in the control group (Figure 4B).

Next, we quantified the bile acid profiles in the liver and feces to evaluate the effect of CUMS on bile acid biosynthesis and metabolism. CUMS significantly increased the secondary bile acid levels in the feces and liver (Figures 3B,C). In particular, the model group liver samples showed increased DCA (148%↑, p = 0.0198), TDCA (166%↑, p = 0.0222), and TLCA (137%↑, p = 0.0028) levels and the model group feces samples showed increased DCA (117%↑, p = 0.0343), TDCA (304%↑, p = 0.0034), and GDCA (306%↑, p = 0.0061) levels.

Since DCA is the microbial metabolic product of TCA, we calculated the relative TCA-to-DCA ratios in the control and model groups to indirectly address the effects of the gut microbiota. The TCA/DCA ratio significantly decreased in the serum and liver of model group (Figure 5). Our results indicated that CUMS markedly promoted intestinal secondary bile acid formation.

We identified a total of 65 bacteria in the intestinal tract from the intestinal content samples at the genus level and summarized the heatmap of the relative abundance in the top 30 genera (Figure 6A). Pearson’s correlation analysis (Figure 6B) indicated that increased secondary bile acid levels in the feces significantly and positively correlated with three members of the phylum Firmicutes: Ruminococcaceae_UCG-010, Ruminococcus, and Clostridia_UCG-014. This result suggested that changes in the secondary bile acid formation might be associated with altered gut microbiota composition in the intestine.