All animal experiments in this study were conducted in the animal facility of the animal center of South China Agricultural University. Specific pathogen-free adult male Sprague–Dawley rats, purchased from Guangdong Medical Experimental Animal Center, were adapted to the laboratory environment for one week (12 h light/dark cycle, 25 ± 1°C, and 55–65% relative humidity with unrestricted access to food and water). The experimental process was presented in Supplementary Figure S1A. After acclimatization, the rats were randomly allocated into two groups: healthy group (HC; N = 12) and chronic unpredictable mild stress group (CUMS; N = 54). The detailed CUMS procedures were listed in our previously published studies (Qu et al., 2019). Briefly, CUMS was induced by eight types of stressors: excessive illumination; cage tilt (45°C, 24 h); wet bedding (24 h); tail suspension (6 min); tail pinch (5 min); food and water deprivation (24 h); bake oven (60°C, 5 min) and flashing light at 150 flash/min for 5 min (Li et al., 2013; Liu et al., 2015). The rats were subjected to randomly select one or two stimuli every day and were not to allowed to repeat the same stimulus for two consecutive days. Each stimulus was guaranteed to be applied to the rats about thirteen times. Overall, the animal experiment lasted 15 weeks. Except for HC rats, all depressed rats received eight different chronic unpredictable mild stimulation for 14 weeks. The body weights of all rats were recorded every week. The CUMS-resistant rats were filtered out by the sucrose preference test (SPT), open field test (OFT), and light/dark test (LDT) at the eight weeks of exposure to CUMS. The remaining rats were administered with Xiaoyaosan (XYS; N = 12), amitriptyline (Ami; N = 6), and fluoxetine (Flu; N = 7) respectively for six weeks by oral administration. Then, after six weeks of treatment with Xiaoyaosan (XYS), amitriptyline (Ami), and fluoxetine (Flu), the depression-like behaviors in rats were evaluated by the SPT, LDT, and FST. Fecal samples from the three individual rats in each group (HC: N = 3; CUMS: N = 3; Ami: N = 3; Flu: N = 3; XYS: N = 3) were selected for metagenomic analysis at 15 weeks and stored at −80°C for subsequent analysis.

DNA was extracted from fecal samples with the QIAmp DNA Stool Mini Kit (Qiagen, USA). Shotgun metagenomic sequencing of all samples was again performed using Illumina HiSeq 4,000 sequencer in Shanghai Majorbio Bio-pharm Technology Co. Ltd. Sequence data associated with this project have been deposited in the NCBI Short Read Archive database (Accession Number: PRJNA910226). Then, the metagenomic raw reads were adjusted quality using Sickle, ¹ and low-quality bases and reads were removed. Burrows–Wheeler Aligner ² compared raw reads with the rat genome to remove the relevant host information. After quality control, the self-written Python script was used for subsequent data processing and analysis. All the clean paired-end reads were de novo whole-genome assembled using SPAdes with the following parameters: -meta -k 21,33,55,77 -only-assembler (Nurk et al., 2017). Scaffolds longer than 1,000 bp were retained for downstream analysis.

To determine confident viral scaffolds, the viral sequence identification standard operating procedure (SOP) with VirSorter2 V.3 from protocols.io ³ was employed, combining VirSorter2, DRAMv, and CheckV, compliance with standard criteria based on viral and host gene counts, score, hallmark gene counts, and contig length (Guo et al., 2021). Additionally, a viral detection pipeline developed by Nayfach et al. (2021) was used to complement the identification of virus sequences. Briefly, this method used a combination of four viral signatures to identify viral genome fragments: (i) the presence of viral protein families, (ii) the absence of microbial protein families, (iii) the presence of viral nucleotide signatures, and (iv) multiple adjacent genes on the same strand. Metagenomic scaffolds were searched for viral and microbial protein families using hmmsearch against Integrated Microbial Genome/Virus (IMG/VR) and Pfam-A database, respectively (Nayfach et al., 2021). VirFinder v.1.1 was applied to evaluate the presence of viral nucleotide signatures. Detailed procedures for selecting viral sequences from the microbiome were performed as previously described (Ren et al., 2017). CD-Hit (Fu et al., 2012) was employed to cluster (95% consistency and 85% coverage) the identified viral sequences, in which the longest gene in each group was selected as the representative sequence, so as to obtain 27,079 DNA non-redundant viral scaffolds (NR_Scaffolds). Because acquiring high-genomes was crucial for our downstream analysis, we retained only NR_Scaffolds with a length of at least 5 kb and defined them as viral operational taxonomic units (vOTUs). CheckV divided the vOTUs into different groups according to assembly quality and completeness for downstream analysis.

The vOTUs were clustered with viral sequences from the IMG/VR v 3.0 database (BLASTn; E value ≤10⁻⁵, nucleotide similarity ≥95%, and covered length ≥ 85%) to be assigned to a known viral taxonomic species. Referring to Nayfach’s research, we clustered again the vOTUs and the IMG/VR database virus sequences using a combination of gene sharing and average amino acid identity (AAI), and determined the taxonomic information of vOTUs at the genus and family levels (Nayfach et al., 2021). Besides, Contig Annotation Tool (CAT) v5.0 (von Meijenfeldt et al., 2019) was also used to taxonomically classify the vOTUs.

To calculate the relative abundance of vOTUs, metagenomic reads of all samples were mapped to the vOTU scaffolds, using Bowtie 2 with default settings (Langmead and Salzberg, 2012). The abundance of each vOTU was calculated by adopting the following equation:

where N_{mapped reads} is the number of mapped reads for each vOTU sequence; L_reads is the sequence length of the metagenomic reads (150 bp in this study); L_{vOTU sequences} is the sequence length of the vOTU; s is the size of each metagenomic data set (Gb).

In specific, the PhaMer software was used to predict the sequences of bacteriophages within the vOTUs identified above (Shang et al., 2022). Then, the phage with a temperate or lytic phage was classified by PhaTYP (Shang et al., 2023). Probability scores attained from PhaTYP were between −1 and 1, which was represent the probability of “Lytic” or “Temperate.”

We used three different methods to predict host-virus connections. (1) Blastn: the vOTUs were aligned with the collected fecal samples metagenomic assembled scaffolds with the following parameters: bitscore ≥50, e-value ≥10⁻³, identity ≥70%, and matching length ≥ 2,500 bp (Emerson et al., 2018), (2) spacer hit: the CRISPR-Cas spacer in the metagenomic sequencing data was predicted using the PILER-CR, ignoring CRISPR-Cas with less than three spacers. BLASTn was used to compare the retain spacers with the vOTUs, and matches were selected when a maximum of one mismatch or gap over ≥95% of the spacer length (Nayfach et al., 2021), and (3) tRNA: Identification of tRNAs from the vOTUs and metagenomic assembled scaffolds with ARAGORN v1.2 (Laslett and Canback, 2004), and compared these tRNAs using BLASTn. Only hits with a 100% identity over 100% of the length were considered (Paez-Espino et al., 2016).

The vOTUs sequence was used as queries for further functional annotation. The open reading frames (ORFs) within these vOTUs were predicted using Prodigal v2.6.3 (Hyatt et al., 2010). Then, the ORFs were used as queries to search for functional genes using BLASTp against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The best hit comparison results of e < 1e-3 were screened to obtain the virus function information.

The co-correlation between the gut bacteria and viruses was constructed by calculating each pairwise Pearson’s rank correlation. Only the bacterial genera and viruses with a strong relation (|ρ| > 0.7 and p < 0.01) were considered. The network layout was visualized using Cytoscape v3. 7. 1 (Shannon et al., 2003).

The behavior results (body weight, SPT, LDT, and FST) were compared between groups using one-way analysis of variance (ANOVA). To evaluate the α-diversity of these fecal samples, we calculated the vOTUs diversity using the Shannon index and richness using the Sobs index. Nonmetric multidimensional scaling analysis (NMDS) based on Bray–Curtis distance was used to assess the variation of viral composition among different groups. Kruskal–Wallis test was performed to compare differences in abundance between each group, which were displayed as means ± SEM. Moreover, the differential KEGG pathway representation was identified using linear discriminant analysis effective size (LEfSe) between the two groups. The LDA score for distinguishing features was set at 2.0. Statistical analyses were conducted using the software SPSS (version 23.0, Armonk, NY, USA), and plots were generated from R and GraphPad Prism (version 8.0).

The experimental scheme was presented in Supplementary Figure S1A, the depression-like behaviors in rats were assessed by the sucrose preference test (SPT), open field test (OFT), and light/dark test (LDT). On week fifteen immediately post-Ami, Flu, and XYS, the body weight of the XYS group was significantly higher than that of CUMS group, whereas no weight gain was significantly observed after Ami and Flu treatment compared to the CUMS rats (Supplementary Figure S1B). After six weeks of treatment with Ami, Flu, or XYS, the Ami, Flu, and XYS rats all demonstrated a significantly increased sucrose preference ratio and rearing and crossing numbers and lowered dark time (Supplementary Figure S1C–E). Among these, the SPT test results of XYS were relatively higher than that of the other two groups, and the Flu group performed best in the OFT experiment and the LDT experiment. These results revealed that antidepressants and XYS could improve depression-like behaviors in CUMS-induced rats.

This study utilized the combination of shotgun metagenomics and virus sequence recognition algorithms analysis to explore the gut virome of the fecal samples from healthy, depressed, and drug-treated rats. Each fecal sample generated about 9 Gb raw reads. After quality control, more than 97% of the raw reads were considered high-quality-filtered reads, so this sequencing dataset contented the requirements of subsequent analysis (Supplementary Table S1). Then the clean reads were assembled into longer scaffolds using SPAdes. SPAdes are well proven to perform well in de novo whole-genome assembled (Nurk et al., 2017). In total, 776,760 scaffolds longer than 1,000 bp were attained with an average length of 3,844 bp. Screening is another crucial process to identify viral sequences from complex scaffolds. We established two bioinformatics methods to screen and identify virus sequences as shown in Supplementary Figure S2. Altogether, in all 5,029 vOTUs were discovered and quantified. We then estimated the quality and length of each viral genome with CheckV (Nayfach et al., 2021). As illustrated in Supplementary Figure S3A, 30 (0.60%) of the vOTUs were classified as complete quality, 378 (7.52%) as high quality, and 4,621 (91.88%) as genome fragments (<90% complete). In addition, nowadays, the largest virus genome found in a human sample is about 393 kb (Zhao et al., 2022). In this study, the genome length corresponding to different completeness was shown in Supplementary Figure S3B, including nine that were ≥ 393 kb (Supplementary Table S2), which may correspond to huge viruses (Al-Shayeb et al., 2020). In particular, one virus (termed vOTU-1) was assembled with lengths of 862 kb, which is longer than the largest phage genome reported so far (735 kb; Al-Shayeb et al., 2020). These results increase our understanding of huge viruses in the gut ecosystem.

To assess the effect of oral antidepressants Ami, Flu, and XYS on gut viral community structure, we calculated each sample’s alpha- and beta- diversity. First, all the identified vOTUs were assessed for alpha diversity using Shannon and Sob indexes. Shannon index showed that only Ami treated mice (Ami: Shannon index = 6.741 ± 0.221 [mean ± SEM]) had reversed viral diversity when compared to CUMS mice (CMUS: Shannon index = 6.553 ± 0.131 [mean ± SEM]), and the Shannon index in Ami group was close to that in HC group (HC: Shannon index = 6.822 ± 0.155 [mean ± SEM]), due to no statistical significance between them (p > 0.05, Figure 1A). Besides, a recovered trend was detected for the viral richness in Ami (Ami: Sobs index = 3962.667 ± 112.036 [mean ± SEM]) and XYS (XYS: Sobs index = 4,024 ± 43.012 [mean ± SEM]) groups compared to that in CUMS (CUMS: Sobs index = 3918.667 ± 151.099 [mean ± SEM]) rats without treatment, but no significant difference was noticed (p > 0.05, Figure 1A).

We then compared β-diversity at the vOTUs level using the Bray–Curtis metric. NMDS demonstrated that the five groups clustered distinctly with a stress value of 0.12 (the value of less than 0.2 indicates the validity of NMDS) (Figure 1B). The NMDS results suggested that the fecal virome of the CUMS rats displayed a distinct composition compared with HC rats (Figure 1B). Intriguingly, with Ami, Flu, and XYS treatment, the viral composition of rats was also separated from CUMS rats, as assessed by NMDS analysis (Figure 1B).

To further quantify the changes in gut virome after treatment with antidepressants and XYS, we attempted to classify vOTUs at the species, genus, and family levels. At the family level, 58.18% of the vOTUs were classified into eight families, which were mainly dominated by bacteriophages from the unclassified_Caudovirales, Siphoviridae, and Myoviridae (Figures 2A,B). Furthermore, our study also found several gut viruses in the giant virus family Mimiviridae (Figure 2B). At the genus level, 75.58% of the vOTUs were unclassified (Supplementary Figure S4A). No vOTUs were annotated at the species level. These results show that there was still a large proportion of unclassified viruses at the family and genus level, suggesting gut unexplored viral diversity.

Next, we performed variation analysis to distinguish differentially present taxa between the five groups at family and genus levels. As shown in Figure 2C, unclassified_Caudovirales was higher in the HC group than in CUMS group (Kruskal−Wallis test, p < 0.05), whereas Podoviridae (p < 0.05) was enriched in the CUMS group relative to the HC group. After six weeks of intervention with three antidepressant drugs, the family unclassified_Caudovirales was found to be increased in the Ami and XYS groups, while only the XYS group exhibited significant differences as compared to CUMS (p < 0.05). We also observed that XYS significantly decreased the abundance of Podoviridae compared to those in CUMS (p < 0.05). Moreover, compared to that of CUMS group, although both oral Ami and XYS could improve the abundance of Siphoviridae, only Ami treatment resulted in statistical significance and achieved levels alike to those of HC (p < 0.05). It is worth noting that the family Myoviridae was significantly underrepresented in the XYS group and reached a significant level with respect to HC (p < 0.05). At the genus level, both Ami, Flu, and XYS upregulated the abundance of g_unclassified_Caudovirales compared with CUMS group (Supplementary Figure S4A). However, only the Ami group reached a remarkable difference (p < 0.05; Supplementary Figure S4B). Compared to CUMS, only Flu treatment resulted in a lower abundance of g_unclassified_Inoviridae (p < 0.05; Supplementary Figure S4B). Ami, Flu, and XYS treatment was also associated with decreased the abundance of Hendrixvirus. Among these, Flu and XYS groups were significantly decreased (p < 0.05; Supplementary Figure S3B). Moreover, the abundance of Phifelvirus was enriched in Ami group relative to CUMS, Flu and XYS groups (p < 0.05; Supplementary Figure S4B). Then in order to identify the common and unique gut virome among the five groups, the gut viral community composition was then compared at the vOTUs level. The Venn diagram displayed that 3,926 vOTUs were shared among the five groups, while 59, 43, 53, 26, and 43 were unique to HC, CUMS, Ami, Flu, and XYS groups (Supplementary Figure S4C).

To outline the effects of antidepressants and XYS on the function of gut virome, we annotated the identified vOTUs set from the gut virome using the KEGG database. In the first-level KEGG pathways, the most abundant functional categories in the different groups were associated with metabolism and genetic information processing (Figure 3A). Among them, the majority of metabolism functions consisted of global and overview maps, carbohydrate and amino acid metabolism (Figure 3B). As for genetic information processing, the proportion of functional genes related to replication and repair was the highest (Figure 3B). Furthermore, we also observed that some genes were directly related to the bacterial hosts, including cellular community: prokaryotes, drug resistance: antimicrobial and infectious disease: bacterial (Figure 3B).

Additionally, to further investigate the function of the drug-related gut viruses, we applied a LEfSe analysis on the third-level KEGG pathway. As illustrated in Figure 3C, 25 different KEGG pathways were identified between HC and CUMS groups (LDA > 2.0). The 9 CUMS-enriched KEGG pathways consisted of amino acid metabolism, carbohydrate metabolism, xenobiotics biodegradation, and metabolism, metabolism of cofactors and vitamins, biosynthesis of other secondary metabolites, replication and repair, aging, and signal transduction. Compared with CUMS group, Ami, Flu, and XYS treatment resulted in 10, 3, and 7 KEGG pathways differences, respectively (Figures 3D–F). Notably, pathway related to xenobiotics biodegradation and metabolism (toluene degradation and fluorobenzoate degradation) was decreased in Ami, Flu, and XYS treatment group, in line with the decreased functions in the HC group (Figures 3D–F). In contrast, pathways of cellular community-prokaryotes (biofilm formation-Escherichia coli) and translation (aminoacyl-tRNA biosynthesis) were enriched in the XYS and Flu groups, separately, in agreement with the enriched functions in the HC group (Figures 3E,F). Moreover, after Ami intervention for six weeks, we also found that pathways associated with translation (aminoacyl-tRNA biosynthesis), carbohydrate metabolism (c5-branched dibasic acid metabolism), and neurodegenerative disease (prion disease) were increased in Ami group, consistent with the increased function in the HC group (Figure 3D).

The interaction between gut viruses and their bacterial hosts dominates the intestinal ecosystem. The gut virome was mainly composed of bacteriophages (phages). Predicting the life cycles of phage could contribute to deciphering their interactions with bacterial hosts. Toward this goal, we assessed the percentage of lytic or temperate phages in the vOTUs. In all, 3,862 phage vOTUs were identified, including 2,555 (66.16%) temperate phages and 1,307 (33.84%) lytic phages (Figure 4A). Notably, temperate phages were the predominant viruses, accounting for more than 65% of the phages’ vOTUs.

Then, we sought to identify the bacterial hosts of the 5,029 vOTUs using three different bioinformatics methods to present how these viruses affect the health of organisms by infecting microbial hosts. First, we testified that 4.2% (212/5,029), 4.0% (202/5,029), 4.5% (224/5,029), 2.4% (122/5,029), and 4.1% (206/5,029) of viruses in the five groups of the gut ecosystem could be connected to a host when using CRISPR spacers strategy, respectively (Figure 4B; Supplementary Table S3). These results indicated considerable CRISPR diversity between bacterial hosts and active infection. Among the spacers in CRISPR arrays record the history of gut viruses infecting bacterial hosts, although most viruses were targeted by the spacer, CRISPR arrays were discovered in only <0.01% of metagenomic scaffolds in each group (Supplementary Table S6), affirming the limited distribution of this anti-viral defense system (Burstein et al., 2016). Then, we were able to assign 381, 397, 417, 341, and 399 viral sequences to bacterial hosts using tRNA in each group (Figure 4B; Supplementary Table S4). Furthermore, to expand the host-virus network, we identified connections between the viral sequences and metagenomic scaffolds in each group based on direct BLASTn matches (bitscore ≥50, e-value ≥10⁻³, identity ≥70%, and matching length ≥ 2,500 bp), resulting in the connection to the bacterial hosts that more than 10% of viral sequences in each group (Figure 4B; Supplementary Table S5). Among them, we were able to identify 123 consistent viral-host links using two or more methods (Supplement Table S7). Overall, these methods identified 25,106 putative host-virus links enabling host assignment to 48.92% (2,460/5,029) of the vOTUs (Supplementary Table S7). Among these vOTUs, a total of 59.15% (1,455/2460) of the vOTUs could be assigned to multiple microbial hosts (Supplementary Table S7). Most predicted hosts were Firmicutes, Proteobacteria, Bacteroidetes, and Actinobacteria (Figure 4C).

In addition, the gut’s relationship between virome and bacteriome is the core of intestine microbiota balance, but the interaction between viruses and bacteria has yet to be well defined. To further link the correlation, the most abundant 50 vOTUs and their bacterial hosts were visualized based on Pearson’s rank correlation (|ρ| > 0.7, p < 0.01) correlations (Figure 4D). Supplementary Table S8 summarizes the detailed co-occurrence between vOTUs and their bacterial hosts. As shown in Figure 4D, the co-occurrence networks revealed that 42 bacterial genera might serve as the potential hosts for the 45 vOTUs. In which 45 vOTUs, 84% (38/45) viruses were linked to more than one bacterial host, indicating that gut viruses also produce cascading effects (Hsu et al., 2019). Of these vOTUs were mainly assigned to unclassified_Caudovirales, Siphoviridae, and Inoviridae families (Supplementary Table S8). Among these 42 genera, 71% (30/42) bacterial genera could be hosts of different viruses, most of which belonged to Firmicutes, Bacteroidota, Proteobacteria, and Actinobacteria phyla (Figure 4D). Interestingly, Muribaculaceae bacterium DSM 108610, belonging to Bacteroidetes, carried more diverse viruses than other genera, indicating that they are the essential microbes dominating the gut ecosystem. In a word, these results suggested that the potential role of these viruses in regulating the bacterial communities in the gut.