The gut virome is associated with stress-induced changes in behaviour and immune responses in mice

Chronic social stress affected behaviour and immune responses in line with previous murine research^26–28. There was a distinct reduction in social interaction behaviour following stress (Fig. 1g). Immobility times in the forced swim test (FST) and, correspondingly, basal corticosterone concentrations were increased in the stress recipient group in addition to a trend of increased corticosterone post FST (Fig. 1h–j). In addition, there was increased inflammatory cytokine release (IL-6 and IL-10) from concanavalin A (ConA)-stimulated splenocytes (Fig. 1k,l), as commonly seen following chronic social stress in mice (Extended Data Table 1a)^29,30. Increased inflammatory cytokine release from splenocyte cultures following stress has been previously linked with anxiety-like behaviour²⁹ and stimulation with ConA indicates humoral-based immune activation capable of driving inflammation.

FVT is a technique where free viral particles, which are dominated by phages, are extracted from the rest of the microbiota and other particulate material. By using phages to target the bacteriome, this treatment both specifically targets the microbiota–gut–brain axis and is safer than comparable faecal microbiota transplants due to the removal of live cellular microorganisms. FVT has previously been used to specifically target the murine bacteriome with phages from another source^11–14. Phage infection of bacterial hosts consists of highly specific targeting in most cases, being governed by strain-level precision. This contributes to the high inter-individual variability of the phageome^22,31–33. We used an autochthonous FVT model since it theoretically provides the highest probability of active phage infections of bacteriome members due to high accordance between phage–host pairs.

We found an assortment of phages present in the FVT preparations, in most cases dominated by members of the class Caudoviricetes and the family Microviridae, although there are large portions of the viral population that are currently unassigned (Fig. 3c). There were 20 differentially abundant bacterial species that were altered by stress but restored when FVT was performed. Fifteen bacterial species were depleted by stress and protected by FVT, while 5 bacterial species were increased by stress and reduced by FVT (Fig. 3d and Extended Data Fig. 3a). Despite differences in individual species in the bacteriome, no changes were found in alpha- and beta-diversity metrics in the bacteriome or virome (Extended Data Fig. 4a–d), indicating that the bacteriome and virome did not undergo a large shift in community profile but rather, in low-abundance members and functionality.

The FVT oral gavage preparation underwent virome sequencing, allowing for the viral inputs to be characterized. Viral host pairs were predicted using a comprehensive collection of viral databases that cover the currently available bacteriophage community; however, these tools are still limited by the vast unknown portion of the bacteriophage virome³¹. Therefore, our bacteriophage sequencing data contain a small percentage of phage contigs that had a predictive host. Predicted hosts of the FVT inocula are reported at the phylum, order and family levels (Extended Data Fig. 3b–d). Phages targeting the phylum Deferribacterota and Bacteroidota (previously Bacteroidetes) were the most common, followed by Bacillota (previously Firmicutes) and to a lesser extent Pseudomonadota (previously Proteobacteria) and Halobacteriota. The most common phages targeted the orders Deferribacterales, Bacteroidales, Lachnospirales, Lactobacillales and Oscillospirales. Lastly, the most common phages targeted the families Mucispirillaceae, Bacteroidaceae, Erypsipelotrichaceae, Lachnospriaceae, Muribaculaceae, OLB10 and Streptococcaceae. However, the majority of the phage–host pairs in the murine microbiome are unknown.

Phage–host pairs from the FVT preparation were then cross-referenced with the differentially abundant taxa (from faecal bacteriome sequencing data) that were affected by stress and restored by FVT. The identified viral taxa that overlapped with the differentially abundant bacterial species were depicted with the predicted phage hosts at the phylum level (Fig. 3e). At the phylum level, 18 of the 20 differentially abundant bacterial species were represented in the predicted phage–host pairs. There were 6 Pseudomonadota, 6 Bacillota, 3 Bacteroidota and 1 each of Verrucomicrobiota, Actinobacteriota and Cyanobacteria pairs. At the order level, there were 2 Burkholderiales, 1 Erysipelotrichales and 1 Bacteroidales pairs. At the family level, there was 1 Lachnospiraceae and 1 Erysipelotrichaeceae pair. This evidence supports bacteriome modulation by FVT, which in turn aided the restorative effects of the phage-derived treatment of stress. Given that phage infections in the gut are overwhelmingly temperate, the phage–host pair associations could be lytic or lysogenic in nature depending on environmental conditions, which can then have direct and/or indirect effects on bacterial host proliferation^9,34. Phage–bacterial cascading interactions are shown as a correlational phage–host network plot (Extended Data Fig. 5). All phage–bacterial correlations depicted show significant association (q < 0.1), and specifically Caudoviricetes 48–Geofilum rubicundum, Unassigned class 38–Kriegella aquimaris and Unassigned class 38–Odoribacter splanchnicus show significant restoration effects following stress and FVT treatment (P < 0.05).

Across the differentially abundant Pseudomonadota taxa, the individual species increased in abundance following FVT compared with stress (Advenella mimigardefordensis, alpha proteobacterium AAP38, Archangium gephyra, Nitrosomonas aestuarii, Paraburkholderia kururiensis and Sinorhizobium medicae), indicating that FVT was a positive regulator. Five Bacillota taxa displayed increased abundance (Bacillus aurantiacus, Bacillus campisalis, Butyrivibrio sp. WCD3002, Cohnella laeviribosi and Erysipelotrichaeceae bacterium NK3D112). Conversely, the Bacteroidota species (Paludibacter sp. 47-17, Pontibacter akesuensis, Thermonema rossianum) decreased in abundance following FVT, as well as a single Bacillota species (Paenibacillus sp. 1-49) and the single Cyanobacteria species (Oscillatoria sp. PCC 10802), indicating that in these cases FVT was a negative regulator. These changes are attributable to stress modulation and either protection or reduction by the FVT treatment. Our data provide evidence that phage–host pairs are present and are potentially involved in the attenuation of stress-coping behaviours within our experimental parameters.

To assess the functional potential of the FVT inocula and bacteriome following the different treatments, gut–brain modules (GBM) and gut–metabolic modules (GMM) were characterized^35,36. In the FVT inocula viral reads, there were three GBMs and ten GMMs identified (Extended Data Fig. 6a,b). Interestingly, the capacity for S-adenosylmethionine synthesis was found in each FVT inoculum and was significantly restored in the bacteriome GBM following stress and FVT. S-adenosylmethionine plays roles in a multitude of biochemical pathways both in the microbiota and across different tissues, including synthesis and metabolism of monoamine neurotransmitters^37–40. We also assessed the restoration of GBM and GMM abundance by FVT in the bacteriome, but no other modules were found to be statistically significant following false discovery rate adjustment (Extended Data Tables 2 and 3).

Finally, to confirm that the FVT could also impact brain circuits underpinning anxiety and stress coping, we performed deep RNA sequencing in the hippocampus and amygdala of all animals. These stress-sensitive brain regions were selected on the basis of previous studies indicating that these regions were susceptible to stress and responsive to modulation by the microbiota^48–50. The top enriched Gene Ontology (GO)⁵¹ terms from genes affected by stress and restored by FVT were selected to specifically target functions directly involved with stress, social behaviour, the virome and the microbiota–gut–brain axis. These GO terms consist of genes associated with fear response, immune system processes, neurotransmitter level modulation, postsynaptic neurotransmitter receptor functionality, stress response, social behaviour and viral processes.

In the first study, 20 male C57Bl/6J mice (n = 10 per group) and 20 male CD1 non-experimental resident aggressor mice were used; in the second study, 30 male C57Bl/6J mice (n = 10 per group) and 50 male CD1 mice were used and housed within the Biological Services Unit, University College Cork. Mice were received at 8 weeks of age from Envigo (UK) and randomly assigned to groups. Animals were habituated for 1 week after which the ‘stress’ groups were single housed while the control groups remained group housed. Food and water were provided ad libitum throughout the study. The holding room was under a 12 h light/dark cycle (light cycle began at 07:00), with a temperature of 21 ± 1 °C and humidity of 55 ± 10%. Body weight was monitored twice per week. To provide an adequate sample size to detect changes; a priori power analysis was performed using the software G*Power (v.3.1) considering a power of 80% (1-beta-error) and a 0.05 alpha-error. Investigators were not blinded during allocation of experiments; however, analyses were conducted in a manner blinded to the experimenter. All procedures were conducted with approval from the Animal Experimentation Ethics Committee (AEEC) at University College Cork and the Health Products Regulatory Authority (HPRA) in accordance with the recommendations of the European Directive 2010/63/EU.

There were two groups, Stress and Ctr, used in the first experiment to test whether chronic social stress changed the bacteriome and virome. Animals (assigned randomly) in the Stress group underwent chronic psychosocial defeat stress which consisted of a semi-randomized sequence of social defeat and overcrowding stressors over a period of 20 d (see Fig. 1a). Behavioural tests were performed between 9:00 and 12:00 as follows: (1) Day 21, social interaction test; (2) Day 22, forced swim test. At the end of the study, animals were euthanized on Day 23 by decapitation and tissues were collected for subsequent analyses.

Animals were randomly assigned to one of three groups used in the second experiment: Ctr, Ctr Stress and FVT Stress. At the start of chronic social stress, animals additionally received treatment by oral gavage; Ctr and Ctr Stress groups received vehicle saline magnesium (SM) buffer (50 mM Tris-HCl, 100 mM NaCl, 8.5 mM MgSO₄; pH 7.5), and the FVT stress group received the purified autochthonous faecal virome from the pre-stress faecal collection. Animals in the ‘Stress’ group underwent chronic social defeat stress 19 times over 22 d (animals were not stressed on the three behaviour testing days). Control animals did not receive stress but were handled evenly by researchers at the same time. Behavioural tests were performed between 9:00 and 12:00 as follows: (1) Day 11, social interaction test; (2) Day 15, elevated plus maze; (3) Day 19, forced swim test (see Fig. 2a). At the conclusion of the study, animals were euthanized on Day 23 by decapitation and tissues were collected for subsequent analyses.

Faecal samples were collected in a PBS-glycerol (20%) solution, frozen (−80 °C) and later used to generate phage-rich material used for FVT oral gavage. Approximately 2.5 g of faecal pellets were resuspended in 10 ml filter sterilized SM buffer and 2 ml 1 M NaHCO₃ (to reduce stomach acidity). The faecal slurry was vigorously vortexed and then centrifuged at 5,000 g for 20 min. The supernatant was then filter sterilized (0.45 μm filter) before undergoing overnight dialysis (SpectraPor Float-A-Lyzer Dialysis Device G2, 100 kDa, Repligen).

In the first experiment, mice underwent 20 d of chronic psychosocial stress, which comprised social defeat stress and/or overcrowding (see Fig. 1a for experimental timeline)²⁶. Experimental mice were stressed over a consecutive 20-d period in an unpredictable social defeat/overcrowding schedule. Social defeat sessions (described below) included exposure to an aggressive CD1 mouse and attack or submissive posture. Overcrowding sessions included co-housing for 24 h in a standard holding cage with n = 10 per cage.

In the second experiment, chronic social defeat stress was carried out daily for 19 consecutive days (excluding behaviour test days; see Fig. 2a for experimental timeline)⁹⁴. Chronic social defeat stress was selected to eliminate potential coprophagy events that may occur during the overcrowding portion of the chronic psychosocial stress paradigm that may confound the effects of FVT treatment. CD1 mice were screened for aggression over 2 d before social defeat sessions. CD1 mice with the shortest attack latencies were selected as the aggressors to be used in subsequent social defeats. For each social defeat session, a CD1 mouse novel to the experimental mouse was used. The session would involve an initial exposure to the CD1 mouse in a clean cage with fresh bedding to reduce coprophagy. The session would last until the first attack, expression of submissive posturing by the experimental mouse, or until 5 min had passed. The latency to attack or submission was recorded. The mice were then separated by a clear Plexiglas perforated wall that allowed non-physical contact for 2 h. After 2 h, the separator was removed, a second defeat session was conducted, then mice were returned to their home cages. Social defeat sessions were carried out in the afternoon during the light cycle. Control mice were handled an equal amount daily in lieu of social defeat.

For study one, animals underwent the social interaction test⁹⁵ after 20 d of chronic psychosocial stress; for study two, animals were tested after 10 consecutive days of chronic social defeat stress. Experimental animals were placed in a rectangular plastic box (31 × 39 cm²) containing an empty wire mesh cage (9.5 × 7.5 cm²) and allowed to freely explore for 2.5 min, then the mouse was removed and placed in its home cage for 1 min while a novel CD1 mouse was placed into the wire mesh cage. The experimental animal was then returned to the social interaction box and allowed to freely explore for 2.5 min while the CD1 mouse was in the wire mesh cage. At the conclusion of the test, the mice were returned to their home cages and the box and wire mesh cage were cleaned with 70% ethanol. Trials were recorded with a ceiling-mounted camera and evaluated using Ethovision 11.5 (Noldus). Social interaction and corner ratios were calculated by taking the ratio of the time spent in the interaction zone (area around the mesh cage) or in the corner (the two furthest points from the mesh cage) during the habituation versus when the CD1 was present.

The elevated plus maze test was used to assess anxiety-like behaviour and was conducted following ref. ²⁸. The maze apparatus was elevated 1 m above ground and consisted of a grey plus (+)-shaped maze with two opposing open arms and closed arms (50 × 5 cm arm length, 15 cm closed wall height, 1 cm open wall height). Experiments and habituation were conducted in red light (5 lux). Mice were habituated to the room for 1 h before the test, then allowed to explore the maze for 5 min. The maze was cleaned with 70% ethanol after each animal trial. Experiments were recorded with a ceiling-mounted camera and videos were scored blind for time spent in open arms and open arm entries (entrance defined as all paws in arm).

The forced swim test was used to assess depressive-like stress-coping behaviour and was conducted according to ref. ²⁸. A glass cylinder (21 cm diameter) was filled with 24 ± 1 °C water to a height of 15 cm; water was emptied and replaced after each individual run. Experimental animals were placed in the water for 6 min while being recorded by a ceiling-mounted camera. After each test, animals were gently dried off with a towel and returned to their home cage. The last 4 min of the 6 min test was used to manually score immobility behaviour from the video recorded by a researcher blinded to the experimental groups. Mice were judged to be immobile when only making movements necessary to keep their head above water.

Plasma from each animal was sampled by a tail tip 5 min before the forced swim test, then 15, 45 and 90 min after the start of the test. The tail was securely restrained and a diagonal incision was made ∼1–2 mm from the tip of the tail. Approximately 40 μl of whole blood was collected in an EDTA-lined capillary tube (Fisher Scientific, 749311) per timepoint. Samples were deposited in an Eppendorf and centrifuged for 10 min at 3,500 g at 4 °C. Plasma was collected and stored at −80 °C for later corticosterone quantification. The remaining blood cells were used for immune flow cytometry analysis.

Samples were analysed in duplicate in a single assay using 10 μl plasma per sample. The threshold of detection was <32 pg ml⁻¹, with coefficient of variation limit of 20% and concentration units expressed in ng ml⁻¹. Light absorbance was read with a multimode plate reader (Synergy HT, BioTek Instruments) at 405 nm.

Flow cytometry was performed on blood samples collected before and 45 min after forced swim test, and samples were processed on the same day. Staining, gating and analysis were performed according to the following methods⁹⁴ and a gating strategy is provided in Extended Data Fig. 7. Blood was resuspended in 10 ml of red blood cell lysis buffer for 3 min (15.5 mM NH₄Cl, 1.2 mM NaHCO₃, 0.01 mM tetrasodium EDTA, diluted in deionized water). Samples were then centrifuged (1,500 g, 5 min) and resuspended in 45 μl staining buffer (autoMACS rinsing solution, Miltenyi 130-091-222, supplemented with MACS BSA stock solution, Miltenyi 130-091-376) for the staining procedure. Then, 5 μl of FcR blocking reagent (Miltenyi, 130-092-575) was added to each sample. Next, samples were incubated with a mix of antibodies at a 2:71 dilution (CD11b+ Miltenyi, 130-113-234; LY6G Miltenyi, 130-117-500; LY6C Miltenyi, 130-111-917; CD62L BioLegend, 104418; MHC-II Biolegend, 107632; CD4 Biolegend, 100548; CD44 BD Biosciences, 563736; CD69 Miltenyi, 130-115-575) and incubated for 30 min on ice. Samples were then fixed in 4% paraformaldehyde (PFA) for 30 min on ice. Fixed samples were resuspended in staining buffer and analysed the following day on the FACSCelesta flow cytometer (BD Biosciences) with FACSDiva (v.8.0, BD Biosciences) and CellQuest Pro (v.5.2, BD Biosciences). Data were analysed using FlowJo (v.10). Cells were gated on the basis of FSC-height and SSC-height; subsequently, LY6C^hi monocytes (CD11b⁺, LY6C^hi), neutrophils (CD11b⁺, LY6C⁺, LY6G⁺) and T helper cells (CD11b⁻, LY6G⁻, CD4⁺) were selected. Next, the proportions of LY6C^hi monocytes and neutrophils that contained CD62L or MHC-II were assessed relative to the respective parent population. Then, the proportions of T helper cells that were either CD69⁺ or CD44^hi were similarly assessed relative to the total parent cell population.

Spleens were collected during tissue collection, trimmed of fat and placed in RPMI medium (R8758, supplemented with 10% fetal bovine serum (F7524, Sigma) and 1% penicillin-streptomycin-glutamine 100x (10378-016, Biosciences)). They were then treated with red blood cell lysis buffer (R7757, Sigma), passed through a 70-μm filter and then resuspended in fresh supplemented RPMI medium. Next, 4 ml of splenocytes were seeded in 6-well plates at a concentration of 16 × 10⁶ cells per well. After a 2.5 h rest period, cells were stimulated with lipopolysaccharide (2 μg ml⁻¹, tlrl-eblps, Invivogen) and concanavalin A (2.5 μg ml⁻¹, C5275, Sigma) for 24 h to induce ceiling antigenicity. Then, supernatant was collected and used for inflammatory cytokine analysis.

Cytokine levels from plasma and splenocyte supernatant samples were collected during euthanasia and quantified using the V-PLEX Proinflammatory Panel 1 mouse kit (MSD, K15048D), measuring IFN-y, IL-6, IL-12/IL-23p40, TNF-a and IP-10 on the MESO QuickPlex SQ 120 SECTOR Imager 2400 (MSD). Cytokine quantification was done according to manufacturer guidelines. Samples below the detection limit of each assay were excluded.

In the first experiment, the QIAamp DNA stool mini kit (Qiagen) was used to extract bacterial DNA using the manufacturer’s handbook (2nd edn 2012) ‘Isolation of DNA from Stool for Pathogen Detection’. Samples were prepared for 16S sequencing using the Nextera XT DNA library prep kit (Illumina), as described in the Illumina 16S library preparation workflow. 16S bacterial rRNA gene was amplified using primers targeting the V3–V4 hypervariable region (Forward: 5′TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG; Reverse: 5′GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC) (Sigma Aldrich). The Illumina V3–V4 primers were selected for their high coverage (94.5% bacteria) while remaining in the amplicon size necessary for sequencing⁹⁶. 16S rRNA amplicons were sequenced on the Illumina MiSeq platform (Teagasc, Moorepark, Ireland).

In the second experiment, The QIAamp Fast DNA stool mini kit (Qiagen) was used to extract DNA from faecal pellets. This was performed following manufacturer’s guidelines. Before extraction, however, 20 μl of the ZymobIOMICs spike-in control 1 (Zymo Research) was added to aid in calculating absolute abundances. DNA concentrations were then measured using the Qubit High Sensitivity assay (Thermo Fisher Scientific) before adjustment to 0.2 ng μl⁻¹ for Nextera XT library preparation (Illumina). Library quality was assessed with the 4200 Agilent Tapestation (Agilent Technologies) using the High Sensitivity D1000 ScreenTape assay (Agilent Technologies) before whole metagenomic sequencing was performed on an Illumina Novaseq platform using 2 × 150 nt paired-end chemistry (Genewiz).

Faecal viromes in the initial experiment were processed using a verified protocol^34,97. Briefly, the samples from control and stress groups were pooled (n = 3–4 per pool) to create 1–2 g faecal samples. Pooled samples were resuspended in 10 ml SM buffer, vortexed for 5 min, cooled for 5 min on ice and centrifuged at 5,000 g for 10 min at 4 °C to remove large particulate and bacterial cells, then supernatants were transferred to new tubes and centrifugation was repeated. Supernatant was filtered twice through a 0.45-µm-pore ceramic filter. Following further refinement and treatment with proteinase K and phage lysis buffer, sample DNA was extracted using the phenol-chloroform method. After DNA extraction, samples were processed for shotgun metagenomic sequencing using the Nextera XT DNA library prep kit (Illumina). Virome metagenomic sequencing was performed on the Illumina MiSeq platform (Teagasc, Moorepark, Ireland).

In the second experiment, an updated protocol described in ref. ⁹⁸ was used. For FVT gavage samples, 400 µl were processed, and to enable comparisons of viral load across samples, 10 µl of 10⁸ plaque-forming units per ml of lactococcal phage Q33 were added to each sample. For faecal samples, a single faecal sample from each mouse was pooled per group in 650 µl of SM buffer. Fresh 0.5 M dithiothreitol was prepared in 1 ml of SM buffer. A volume of 16 µl of the dithiothreitol stock was added to samples to achieve a final concentration of 12 mM and samples were incubated at 37 °C for 30 min. Centrifugation at 4,000 g for 30 min at room temperature was used to pellet debris and bacterial cells. Subsequently, 400 µl of liquid was removed and treated with 40 µl of DNase/RNase buffer (50 mM MgCl₂, 10 mM CaCl₂), 12 µl (24 U) of Turbo DNase (Thermo Fisher) and 4 µl (40 U) of Fermentas RNase I (Fisher Scientific) at 37 °C for 1 h. Nucleases were inactivated by incubating at 65 °C for 10 min. With free nucleic acids now removed, virus-like particles (VLPs) were lysed using the QIAgen blood and tissue kit following manufacturer guidelines. However, an elution volume of 30 µl of Qiagen AE elution buffer (10 mM Tris-HCl, 0.5 mM EDTA, pH 9.0) was used to increase the final concentration of nucleic acid obtained.

A reverse transcription reaction was performed using the SuperScript IV First Strand Synthesis System (Invitrogen/Thermo Fisher Scientific) with 11 µl of purified VLP nucleic acid sample and random hexamer oligonucleotides according to manufacturer protocol. The concentration of DNA purified using a DNeasy blood and tissue kit (Qiagen) was determined using the Qubit dsDNA HS kit and the Qubit 3 fluorometer (Invitrogen/Thermo Fisher Scientific). DNA/cDNA yields varied between 0.05 and 29 ng µl⁻¹, with some samples being below the detection limit. Before library preparation, 22 µl of each sample underwent reverse transcription using the Superscript IV kit (Invitrogen/Thermo Fisher Scientific) according to manufacturer protocol to include any RNA viruses in the analysis.

Subsequently, samples were sonicated following adjustment of the volume to 52.5 µl with low-EDTA TE buffer. Shearing of unamplified DNA/cDNA mixture was performed on an M220 Focused-Ultrasonicator (Covaris) with the following settings: peak power of 50 W, duty factor of 20%, 200 cycles per burst, total duration of 35 s. Samples were then concentrated using the Genomic DNA Clean and Concentrator-10 kit (Zymo Research) and eluted in 16 µl aliquots. Library preparation was then carried out using the Accel-NGS 1S Plus kit (Swift Biosciences) according to manufacturer instructions. However, a final and additional bead clean-up step of 0.8 DNA/AMPure beads v/v ratio was included. Library QC to determine fragment length distribution and quantitation was performed with the 4200 Agilent Tapestation (Agilent Technologies) using the High Sensitivity D1000 ScreenTape assay (Agilent Technologies) and Invitrogen Qubit. The normalized pooled library was then sequenced using 2 × 150 nt paired-end sequencing chemistry on the Illumina HiSeq platform at Genewiz.

We performed quality checks on raw sequences from all faecal samples using the FastQC programme with a quality score of 30 as a threshold. Shotgun metagenomic sequencing data were then cleaned and host genome sequences were filtered using Bowtie2 via the Kneaddata wrapper programme with following parameters: ILLUMINACLIP:/NexteraPE-PE.fa:2:30:10, SLIDINGWINDOW:5:25, MINLEN:60, LEADING:3, TRAILING:3. Reads were aligned to the Web of Life database using Bowtie2, and taxonomic and functional profiling of the microbial community was performed using Woltka. Next, the uniref90-based gene abundance matrix was further collapsed using KEGG Orthology (KO) term mapping via the ‘woltka tools collapse’ function provided within Woltka. Woltka SOP is available online (https://github.com/qiyunzhu/woltka/blob/master/doc/wolsop.sh↗). GBMs³⁵ and GMMs³⁶ were calculated using the R version of the Gomixer tool.

Virome reads were processed to remove adaptors, trim low-quality bases, and filter reads for quality and length using fastp (v.0.23.2)⁹⁹ with the flags ‘-length_required 60–detect_adapter_for_pe –cut_front –cut_tail –cut_window_size 4 –cut_mean_quality 20–trim_front1 10’. Eukaryotic reads were then identified and removed by mapping reads to the Mus musculus and Homo sapiens genomes (GCF_000001635.27_GRCm39 and GCF_009914755.1_T2T-CHM13v2.0, respectively) using Kraken (v.2.1.2)¹⁰⁰. At each step, read quality was assessed using FastQC (v.0.11.9)¹⁰¹.

Reads passing the trimming and filtering steps were then assembled on a per-sample basis using MEGAHIT (v.1.2.9)¹⁰² with default parameters. The resulting contigs from all samples were then pooled and filtered to retain only those contigs ≥3 Kbp in length. Redundant contigs were then removed following an all-vs-all BLASTn (v.2.12.0+) comparison (e-value cut-off of ≤1 × 10⁻²⁰) where the shortest contig was removed for all pairs of a hit for contigs that shared ≥90% identity over ≥90% of the shortest length of the contig.

Putative viral contigs from the non-redundant contig pool were then identified using the following criteria: first, contigs were compared to the viral nucleotide sequences within the National Center for Biotechnology Information (NCBI) RefSeq database release 211, the Gut Virome Database¹⁰³, the Joint Genome Institute (JGI) IMG/VR v.3 (release #5; 12-10-2020)¹⁰⁴ and our in-house crAss-like bacteriophage genome database using BLASTn (v.2.12.0+) with an e-value cut-off of ≤1 × 10⁻¹⁰. For these nucleotide comparisons, any contig found to have ≥85% of its length covered by any combination of hits with at least 50% identity to sequences within these databases was deemed as a potential viral sequence. As phages are particularly heterogeneously distributed among microbiomes if a feature is absent (or unobserved) in approximately half of samples, inclusion in subsequent analysis greatly increases the risk of spurious associations¹⁰⁵. Potential viral contigs were also identified on the basis of protein homology. For this, open reading frames (ORFs) were called within each contig with Prodigal (v.2.6.3)¹⁰⁶ using the ‘-meta’ option, with the resulting ORFs compared to the NCBI RefSeq viral protein database (release 211) using BLASTp (v.2.12.0+; e-value cut-off of ≤1 × 10⁻¹⁰). ORFs were also compared to the hidden Markov models (HMMs) of the Prokaryotic virus Remote Homologous Groups (PHROGs; release v.3) database¹⁰⁷ and the JGI’s inovirus protein families (iPFs) database¹⁰⁸ using hmmscan (HMMER v.3.3.2; e-value cut-off of ≤1 × 10⁻⁵). Contigs were denoted as potentially viral from the results of protein comparisons if the contig encodes for at least 2 ORFs that were hit by any combination of the protein comparisons where this accounted for at least 50% of the total of the contig ORF count. Two stand-alone viral identification tools were also used to identify putative viral contigs from the contig pool, namely VirSorter2 (v.2.2.3)¹⁰⁹ (with ‘–include-groupsdsDNAphage,RNA,ssDNA’) and DeepVirFinder (v.1.0)¹¹⁰ (with default paramaters). Finally, any circular contigs identified by the ViRal and Circular content from metAgenomes tool (VRCA; https://github.com/alexcritschristoph/VRCA↗) were also included in the list of potentially viral contigs.

Contigs that met any of these viral identification criteria were then examined for the presence of rRNA genes using barrnap v.0.9 (https://github.com/tseemann/barrnap↗) with default parameters, with any contigs found to contain rRNA hits then removed from the list of potential viral contigs. Finally, the potential viral contigs were then examined for quality using CheckV v.1.0.1 and its sequence database (v.1.2)¹¹¹. In cases where a contig was identified as proviral by CheckV, the contig sequence was refined to that within the viral start and end positions suggested by CheckV for that contig. As this refinement altered some of the viral contig lengths, contigs where again re-filtered to remove contigs ≤3 Kbp in length.

From the list of viral contigs, the assembled contig corresponding to the lactococcal phage Q33 used as an internal viral standard was identified by BLASTn of the phage Q33 genome (GenBank accession: JX564242.1↗) and the corresponding contig replaced with the reference genome sequence for Q33. Viral contigs were then clustered into viral clusters using vCon- TACT2 (v.0.11.1)¹¹² with the options ‘–rel-mode Diamond –db ProkaryoticViralRefSeq207-Merged–pcs-mode MCL –vcs-mode ClusterONE’. Viral family information was also assigned to contigs using Demovir (https://github.com/feargalr/Demovir↗). Where vContact2 was able to cluster viral contigs with a known virus, the taxa information of the known virus was assigned to the contig, otherwise the information from Demovir was used. The lifestyle of each viral contig (that is, virulent vs temperate) was predicted using BACPHLIP (v.0.9.6)¹¹³. If BACPHLIP assigned a lifestyle with a confidence value of ≥0.95, that lifestyle was assigned to the contig, otherwise, the lifestyle of the contig was denoted as ‘unassigned’. Putative bacterial host assignments for each contig were achieved with iPHoP (v.1.1.0)⁷² using the ‘Sept_2021_pub’ host database and default ‘predict’ parameters.

Quantitative analysis of the faecal virome inocula was achieved by first mapping the per-sample reads passing trimming and filtering onto the final database of putative viral contigs using Bowtie (v.2.4.5)¹¹⁴ in ‘–end-to-end’ mode. Read coverage and mapping counts were then collected using SAMTools (v.1.15.1)¹¹⁵. Spurious read mapping counts to viral contigs were eliminated from the count data by setting mapped read counts to 0 for any viral contig found to have 1x coverage over <75% of its length. Putative viral contigs were included in the quantitative analysis only if: the contig had been determined to have ≥50% completeness by CheckV, the number of CheckV-identified viral genes was ≥1 and the number of genes classified as ‘viral’ by CheckV was ≥10% of the total gene count. For all samples, the abundance of each viral contig in a sample was estimated by normalizing the number of reads mapping to a contig by the length of the contig, followed by normalization for library size via ‘total-sum scaling’¹¹⁶ (that is, length-normalized read count / total number of reads used for mapping in a sample). Absolute viral counts were calculated by comparing the relative abundances of normalized read counts to that of the spiked lactococcal phage Q33 contig. GBMs and GMMs were screened for within the assembled virome sequences by first annotating viral genes using eggNOG-mapper v.2.1.11, then searching the output KO against the GBM database.

Quantitative analysis of the total viral sequence information within the bulk metagenomic sequencing data followed a similar process used in the analysis of viral gavage samples. Namely, reads that passed trimming and filtering by fastp for each sample were mapped to the database of putative viral contigs obtained from the gavage samples as described above. Mapping results were then converted to total-sum scaled read counts as described above. Again, only those viral contigs that passed the CheckV filtering conditions described earlier were included in read count table.

Further data handling of microbiome measures was undertaken in R (v.4.2.2) using the Rstudio GUI (v.2022.7.2.576).

Interplate biases were recalibrated using the Metacal package in R¹¹⁷. Therefore, taxa with a minimal count of <5 were excluded from the bacterial analysis and taxa with a prevalence of <50% were excluded from the phage analysis, as ratios are invariant to subsetting and this study employs compositional data analysis techniques^118,119. This approach was used since phages are heterogeneously distributed among microbiomes and if a feature is absent (or unobserved) in approximately half of samples, then inclusion in subsequent analysis greatly increases the risk of spurious associations¹⁰⁵. Principal component analysis was performed on centred log-ratio transformed (clr) values. Zeroes were replaced using the ‘const’ approach described in ref. ¹²⁰. Beta-diversity was computed in terms of Aitchison distance or Euclidean distance between clr-transformed data and assessed by permutational multivariate analysis of variance (PERMANOVA) using the vegan package. Alpha-diversity was computed using the iNEXT library¹²¹.

Differential microbial abundance was assessed using generalized linear mixed-effects models using the lme4 framework¹²². Custom functions and scripts are available online at https://github.com/thomazbastiaanssen/Tjazi↗. Where appropriate, post hoc pairwise comparisons were performed using one-tailed Tukey’s procedures with a P value cut-off of 0.05, in addition to a P value for the orthogonal restoration comparisons with a threshold of 0.05 and a Benjamini–Hochberg adjusted q-value cut-off of 0.2.

Restoration of microbial taxa was assessed using orthogonal contrasts in a generalized linear mixed-effects model, weighing the ‘Ctr’ and ‘FVT Stress’ group against the ‘Ctr Stress’ group. Figures were generated using ggplot2, ggforce and patchwork in R.

VLPs in the FVT solution were stained with SYBR gold as adapted from a previous method¹²³. Briefly, the FVT solution was diluted and passed through an all-glass vacuum filter onto an Anodisc (25 mm diameter, 0.02 μm pore size; Whatman) for 60 s. Then, the Anodisc was stained with SYBR gold for 20 min, antifade was added and it was then allowed to dry in the dark. Images were taken on a fluorescence upright microscope (Olympus BX53 Upright Research Microscope) at ×100 magnification, and VLPs in images were quantified using a Python script (available at https://zenodo.org/records/10460215↗).

Hippocampus and amygdala RNA was extracted using the mirVANA kit (Qiagen). Messenger RNA sequencing was conducted by Source Bioscience on the Illumina NovaSeq 6000 (S4 flow cell, 150 bp paired-end lane yielding 2.5 billion read pairs). A reference genome of obtained sequences was created using the reference annotation: Mus musculus (organism), GRCm38, UCSC, genome browser (GRCm38.p6), Ensemble¹²⁴. Reads were assessed and filtered for quality with fastqc using default parameters. Genes were annotated and counted with kallisto¹²⁵ using default parameters.

Further handling of transcriptomic data was undertaken in R (v.4.2.2) using the Rstudio GUI (v.2022.7.2.576). Principal component analysis was performed on clr-transformed values. Zeroes were replaced using the ‘const’ approach described in ref. ¹²⁰. Differences in overall gene expression were computed in terms of Aitchison distance (Euclidean distance between clr-transformed data) and assessed by PERMANOVA using the vegan package.

Differential gene expression was assessed using generalized linear models employing the base R stats lm framework. Custom functions and scripts are available online at https://github.com/thomazbastiaanssen/Tjazi↗ refs. ^126,127. Before statistical testing, hippocampal gene transcripts were filtered by variance using the genefilter library in R¹²⁸, with 0.5 set as the cut-off. Restoration of gene expression was assessed using orthogonal contrasts in a generalized linear model, weighing the ‘Ctr’ and ‘FVT Stress’ groups against the ‘Ctr Stress’ group for the hippocampus and amygdala individually.

To assess whether a subset of genes were altered by stress and restored to control levels post FVT, we employed the enrichR framework¹²⁹ on the GO database¹³⁰. In parallel, we identified seven terms in the gene ontology database that corresponded with the objectives of the study design. These seven terms were additionally assessed for enrichment using the base R stats phyper implementation of the hypergeometric test. Figures were generated using ggplot2 and patchwork in R.

16S rRNA, metagenomic and transcriptomic data were analysed using R, and the remaining statistical analyses were carried out using SPSS v.28. Measurements were taken from distinct samples and not measured repeatedly. In outputs where multiple technical replicates were performed (for example, corticosterone and inflammatory cytokine analysis), the mean per individual was used. The restoration by FVT of features altered by stress was analysed by planned orthogonal contrast (P), followed by Benjamini–Hochberg’s adjusted false discovery rate (BH-FDR) of the calculated P value (q) and Tukey’s post hoc test (P; pairwise comparison of Ctr–Ctr Stress, followed by comparison of Ctr Stress–FVT Stress). For the phage–host network (Extended Data Fig. 5), the Anansi package¹³¹ was used whereby data were compared using linear models, followed by BH-FDR and Pearson’s correlation coefficient to plot the phage–host pairs. Unless otherwise noted, significant effects are denoted by: *P < 0.05, **P < 0.01, ***P < 0.001 and †q < 0.2. Technical outliers were removed before statistical analyses (for example, animals that did not complete tests, samples not available for assays, values outside limits of detection and so on). Values that were 2.5 standard deviations from the mean were considered statistical outliers and excluded from the related analyses (maximum of 2 statistical outliers per test). Data distribution was assumed to be normal, but this was not formally tested. Analyses were conducted in a manner blinded to the experimenter.

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary

Source data Supplementary Source Data Supplementary source data.