The bacterial cell lines Lactobacillus plantarum ATCC 793 (Lp793) and Bifidobacterium longum ATCC 15707 (Bl15707) were cultivated from frozen stock in Man-Rogosa-Sharpe (MRS) media and MRS with 0.05% cysteine, respectively, in an anaerobic incubator at 37°C as previously described (28).

C57BL/6J male mice were purchased from The Jackson Laboratory (age 8 weeks; Bar Harbor, ME, USA) and group housed in the centralized animal care facility of the Center for Comparative Medicine and Surgery at the Icahn School of Medicine at Mount Sinai. Male mice were exclusively used in this study to normalize the behavioral effects due to stress; further studies will be conducted to directly compare sex in stress-induced behavioral deficits. All animals were maintained on a 12-h light/dark cycle in a temperature-controlled (20 ± 2°C) vivarium with access to food and water ad libitum. All procedures, protocols and behavioral experiments were approved by the Mount Sinai Institutional Animal Care and Use Committee (IACUC).

Animals were fed a polyphenol-free diet for the duration of all experiments ( Table S1 ). Following the 2-week stabilization and acclimatization period for the mice, animals were placed on their respective treatment for 2 weeks prior to starting the stress protocol. The Bioactive Dietary Polyphenol Preparation (BDPP) was comprised of 1% w/v grape seed polyphenol extract (GSPE; Healthy Origins), 1% w/v resveratrol (BulkSupplements.com) and a 5% w/v concord grape extract (AA Pharmachem, San Diego, USA) made in sterile water. All tested compounds were analyzed by liquid chromatography–mass spectrometry and archived as previously reported (28, 31) in compliance with National Institutes of Health, National Center for Complementary and Integrative Health (Bethesda, MD, USA) product integrity guidelines. Probiotic doses for mice were prepared by growing cultures in bulk as previously described (28). The bacteria were incorporated into the animals’ drinking water at a final dosage of 1.0x10 (9) CFU/day per bacterium, calculated based on the average daily water consumption per cage, which was recalculated and replaced daily. The synbiotic was composed of BDPP and the probiotics and was replaced daily. Bacteria cultures were prepared in bulk, aliquoted at a known density and frozen at -80°C until needed. This viability of this preparation was tested by freezing bacterial stocks as indicated and replating from thawed samples and preforming colony counts.

Mice were randomly subdivided into 8 treatment groups including untreated control (Control), BDPP-only (BDPP), probiotic or synbiotic and each treatment group was treated with one condition, either stressed or non-stressed. During the stress protocol, mice were submitted to 28 days of random mild stressors (CUS), allowed to rest for 28 days (CUS+Rest) followed by re-stimulation to a subthreshold mild unpredictable stress (US) for 7 days (CUS+US) to model recurrent stress. The CUS+US timepoint was controlled to a 7 day US without prior stress exposure as demonstrated in previous studied ( Figure 1A ) (15). For each treatment, condition (stressed or non-stressed) and timepoint, the group size included n=16 animals, combined over three independent experiments where tissues were allocated towards the biochemical and imaging studies equally. This group size was determined based on power analyses from previous data conducted by us (28) and others (15). During the stress periods, animals in all study groups, including the mice for the immunophenotyping, were subjected twice daily to random mild stressors ( Table S2 ). Stressors included 45° cage tilt for 12 h, wet bedding for 10-12 h, no bedding for 10-12 h, food and/or water deprivation for 12 h, 4°C cold exposure for 1 h, cold water swim for 5 min, cage shaking for 20 min, reversed light schedule, restraint stress for 1 h, predator scent exposure for 8 h, or crowding with 12 animals/cage for 1 h. Consistent with our ongoing and published (28) experiments, no significant changes between groups were observed for weight, water or food consumption throughout the testing period. Behavioral assessment on all animals (n =16 per group) was conducted following each timepoint. From the 16 animals per group, 6 mice were perfused and used for immunofluorescent analysis, 6 mice were used for tissue RNA, protein (brain and peripheral tissues), blood and feces collection and the remaining 4 for other purposes not relevant to this manuscript. Animals were sacrificed immediately following behavior, blood was collected by cardiac puncture into heparinized tubes while tissues were frozen on dry ice and stored at -80°C until analysis.

Behavioral experiments were performed with a NIR camera and measured with ANY-maze™ tracking software (Stoelting Co., IL, USA. Version 5.1 Beta). All animals were handled for 5 min per day for 3 days prior to behavioral testing and were habituated to the testing room for 1 h at the beginning of the test day. Details for anxiety-like behavior with the open field test and depressive-like behavior with the forced swim test can be found as previously described (28).

Fecal pellets from six mice were isolated from the distal colon from 6 mice per group. The Qiagen PowerFecal Pro DNA kit (Cat No. 12830-50) was used for DNA extraction with the Omni International BeadRupter 24 bead mill homogenizer. DNA from extracted samples was amplified using Invitrogen’s AccuPrime High Fidelity kit using primers containing adapters for MiSeq sequencing and single-end barcodes allowing pooling and direct sequencing of PCR products (515F/806R) (33). The 16S rRNA gene sequencing was conducted by Diversigen (Houston, TX) (34). The 16S rRNA V4 gene region was amplified from the extracted community DNA by PCR and sequenced with the MiSeq platform (Illumina) using the 2x250 bp paired-end protocol. This generated paired-end reads that overlapped almost completely. Sequences that passed chimera slaying post sequencing were merged and blasted against the 16S specific curated Silva database (v132).

Previous studies showed promising changes in immunological gene regulation in the ileum, spleen and multiple brain regions (28). To confirm these changes in the current study for multiple time points, RNA from the spleen and ileal tissues were extracted using the Trizol reagent (ThermoFisher) while RNA from brain regions was extracted using the RNeasy Mini kit (Qiagen), with the same n = 6 animals per group as used for the 16S sequencing. cDNA synthesis was conducted with the High-Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (ThermoFisher). qPCR was conducted in collaboration with the Quantitative PCR CoRE at the Icahn School of Medicine at Mount Sinai using an ABI 7900HT Real-Time instrument and SDS software. Relative gene expression was assessed using the delta-delta CT method (35) and all primer pairs and annealing temperatures are provided in Table S3 . All data was expressed in relative terms to the housekeeping gene gapdh following its demonstration of stability across tissue and treatment (36).

Protein from brain regions, spleen and ileum from the same mice used for RNA and fecal collection (n = 6) were homogenized in RIPA buffer (Sigma) with an added protease inhibitor cocktail (Sigma) and phenylmethylsulfonyl fluoride (PMSF) phosphatase inhibitor (ThermoFisher). ELISA assays for the cytokines IL-1β (limit of detection, LOD 15.6 pg/ml, R&D Biosystems), IL-10 (31.2 LOD, R&D Biosystems), IL-6 (LOD 4.0 pg/ml ThermoFisher) and IL-17A (ThermoFisher, LOD 4 pg/ml) were conducted as per company instructions.

At the time of sacrifice, mice (n = 6 per group, separate from the RNA, protein and fecal analysis) were cardiac perfused first with saline and then with 4% paraformaldehyde (PFA) and the brain including the brain stem was removed. Brain tissue was preserved in 4% PFA overnight and subsequently transferred to a saline solution containing sodium azide (0.02%). Sagittal sections were cut to 50 µm thickness using a vibratome (Leica VT1000S) and stored in saline with sodium azide. For immunostaining, slices were washed with PBST (PBS + 0.1% Triton X-100) and blocked in 5% normal goal serum. Slices were incubated with primary antibody (rabbit anti-mouse Iba1, clone EPR16588, abcam and rat anti-mouse CD68, clone FA-11, BioRad) and then secondary antibody (goat anti-rabbit AF568 and goat anti-rat AF488, ThermoFisher) before mounting with ProLong Diamond Antifade mounting media (ThermoFisher). Images were acquired on a Zeiss LSM880 Airyscan confocal microscope under an X20/0.8 NA air immersion 20 X objective. Image analysis was conducted with Image J. Microscopy image analysis was performed at the Microscopy CoRE at the Icahn School of Medicine at Mount Sinai.

Samples processed for CyTOF were derived from a different cohort of mice than the reported behavioral and biochemical analyses with n = 3 animals per group. Blood was drawn by cardiac puncture and placed in EDTA blood collection tubes on ice. PBMCs were extracted by serially lysing red blood cells (RBC) in RBC lysis buffer. The spleen and ileum (distal 5 cm of the small intestine) were removed and placed in ice-cold RPMI supplemented with FBS (2%) and HEPES (15 mM). The spleen was washed, macerated and placed in a digestion buffer containing RPMI supplemented with FBS (5%), DNase I (0.5 mg/mL) and collagenase IV (0.4 mg/mL) at 37°C for 30 min with agitation. After 30 min, the spleen was pulverized with a 18G needle and filtered through 100 µM mesh on ice. Ileum samples, with Peyer’s Patches removed, were washed in HBSS to remove fecal matter and cut longitudinally so mucus and remaining fecal matter could be removed. The whole ileum was placed in a dissociation buffer of RPMI supplemented with FBS (5%), EDTA (5 mM) and HEPES (15mM) for 20 min at 37°C with agitation. Following the 20 min, the dissociation buffer and ileum was filtered through a 100 µM filter the lamina propria was transferred to a digestion buffer of RPMI supplemented with FBS (5%), DNase I (0.5 mg/mL) and collagenase VIII (0.4 mg/mL) and incubated 25 min at 37°C with agitation. Like the spleen, the ileum was pulverized with a 18G needle and filtered through a 100 µM mesh on ice.

Samples were processed by the Mount Sinai Human Immune Monitoring Center. Cell counts were performed on the Nexcelom Cellaca Automated Cell Counter (Nexcelom Biosciences, Lawrence, MA, USA) and cell viability was measured utilizing Acridine Orange/Propidium Iodide viability staining reagent (Nexcelom). Cell counts were normalized such that four million cells were taken for downstream processing. After washing cells, live-cell barcoding was performed. Live-cell barcoding allows for the pooling of samples prior to CyTOF antibody labeling, thus eliminating batch staining variability between replicates (37). In this study, cadmium conjugated CD45 and MHCI (H-2) antibodies were utilized for live-cell barcoding, and replicates of each condition were live-cell barcoded together. Fc receptor blocking (Biolegend Inc., San Diego, CA, USA) and Rhodium-103 viability staining (Fluidigm) were performed simultaneously with live-cell barcoding. Next, surface staining (antibody list Table S4 ) was performed by resuspending each pooled sample in a scaled amount of CyTOF antibody cocktail. Next, the eBioscience Foxp3/Transcription Factor Staining Buffer Set (Thermo Fisher Scientific, Waltham, MA, USA) was used to perform intranuclear staining. Samples were then washed and palladium barcoding of each condition was performed utilizing the Fluidigm Cell-ID 20-Plex Pd Barcoding Kit (Fluidigm) following manufacturer’s instructions. Conditions with unique combinatorial palladium labels were pooled such that only three samples remained (PBMC pool, spleen pool, ileum pool). This second tier of barcoding allowed intracellular staining to be performed on all of the conditions/replicates per tissue type in bulk. Heparin blocking (100 units/mL) was utilized to prevent non-specific binding of intracellular antibodies to eosinophils (38). Samples were fixed in 2.4% PFA. 125nM Iridium-193 (Fluidigm) and 2nM Osmium tetroxide (EMS) cell labeling was performed simultaneously with sample fixation (39). Samples were washed twice with CSB and stored in FBS + 10% DMSO at -80°C until acquisition (40).

Prior to data acquisition, samples were washed and resuspended at a concentration of 1 million cells per ml in Cell Acquisition Solution containing a 1:20 dilution of EQ Normalization beads (Fluidigm). The samples were acquired on a Helios Mass Cytometer equipped with a wide bore sample injector at an event rate of <400 events per second (41). After acquisition, repeat acquisitions of the same sample were concatenated and normalized using the Fluidigm software, and palladium-based debarcoding was performed utilizing the CyTOF debarcoding software made available from the Eli Zunder lab (42). These debarcoded files were uploaded to Cytobank for manual data clean-up/live-cell debarcoding.

For manual cleaning, immune cells were first identified based on Ir-193 DNA intensity and CD45 expression; Ce140+ normalization beads, CD45-low/Ir-193-low debris and cross-sample and Gaussian ion-cloud multiplets were excluded from downstream analysis. After this data cleanup, manual gating was utilized to debarcode the live-cell cadmium barcoded replicates. After this final debarcoding, the files were split by population such that each FCS file contained fully cleaned and debarcoded data from its corresponding sample. Lastly, for unbiased analysis of cell-subsets within each sample, the debarcoded files were run through the Astrolabe Diagnostics Data Processing pipeline (43). viSNE plots and cell type gating were preformed using Cytobank (44). viSNE clustering was performed on 22 parameters ( Table S5 ) where equal event sampling was selected using 6666 events in the spleen, ileum and PBMCs, the lowest common denominator in all samples.

All linear regression and Multivariate Adaptive Regression Splines (MARS) analyses were performed in R (version 4.0.2). The samples included were all derived from biological material from the same mice used for the 16S sequencing, RNA, protein and behavioral tests with n = 6 per group. Basic linear regression was performed in R using the lm function (“stats” R package, version 4.0.2). The association of individual bacterial species, gene expression and cytokines with subject behavior was examined with the following formula: log10(response) ~ Behavior + Day + Condition, where Day and Condition served as control variables. MARS analysis was performed in R using the “earth” package (version 5.1.2). The MARS approach allowed for the systematic identification of pertinent primary, secondary or higher levels of interactions among predictor variables. Briefly, MARS constructs a linear model with all possible basis terms (defined as hinge functions of predictors and the products of them), then by performing a stepwise term deletion in the full model, identifies the set of model terms that yield the best Generalized Coefficient of Variation (GCV). To assess the influence of microbiome on behavior, the MARS algorithm was run with the behavior as the response variable, while all available microbiome genera and conditions (i.e. time, stress or treatment) as predictors. To assess the influence of measured genes and cytokines on behavior, the MARS algorithm was run with the behavior as the response variable, while all available genes, cytokines, and conditions (i.e. time, stress or treatment) were predictors. The optimal parameters for MARS analysis, including the maximum degree of interactions, was determined using the generalized R-squared metric, which is the estimated model performance on unseen datasets. To systematically dissect which minor cell population was varying with respect to treatment and/or time, the MARS algorithm was repeated with the cell frequencies jointly as the predictor variable, and treating the timepoint, treatment and stress conditions as the response variable.

All statistical analyses were completed in Graphpad Prism version 8.0. All within group comparisons for the behavior, gene and cytokine data were conducted with 2-way ANOVAs with Tukey’s posthoc test while one-way ANOVAs with Tukey’s post-hoc analysis for the microglia activation marks and immune cell frequency data. For the gut microbiome analysis, between group statistics of the phyla and genera were conducted with the Mann-Whitney U test (due to the non-parametric nature of gut microbiome data) with false discovery rate (FDR) corrections using the Benjamini-Hochberg method. Alpha diversity, representing microbiome variance within a sample, was calculated with the Fisher diversity index and Observed OTU abundances. Beta diversity, between sample variance, was estimated quantitatively with weighted UniFrac analysis distance matrices. Variation in community structure was assessed with permutational multivariate analyses of variance (PERMANOVA) between groups with visualization of the data using Principal Coordinates Analysis (PCoA).

Chronic and recurrent stress were modelled using the CUS paradigm ( Figure 1A ). This paradigm was used to elucidate the gut microbiome-peripheral-brain-inflammation connection in driving chronic- and recurrent-stress induced behaviors. As previously shown, behavioral changes were not observed following US (15); however, significant depressive- ( Figure 1B ) and anxiety-like ( Figure 1C ) behaviors were observed in the stressed vehicle controls following CUS. Both behavioral impairments recovered following the rest period (CUS+Rest) and were re-observed with CUS+US re-stimulation. The probiotic and synbiotic attenuated depressive-like behavior compared to stressed vehicle controls following CUS, while the synbiotic rescued the phenotype following CUS and CUS+US ( Figure 1B ). All treatment groups elicited a beneficial effect on anxiety-like behavior following CUS and CUS+US compared to the stressed vehicle control ( Figure 1C ). Notably, only BDPP and the synbiotic, however, improved anxiety-like behavior compared to their respective unstressed controls at the CUS timepoint. Taken together, the synbiotic promoted more consistent and robust resilience to stress-induced behavioral impairments compared to its components.

Studies have shown that microglia activation in limbic brain regions mirror behavioral responses of chronic and recurrent stress (15). In this study, microglia activation was determined by the activated CD68 positive surface area. There were no variations in microglia activation or number in the prefrontal cortex (PFC) with US, but following CUS, microglial activation was elevated in the stressed vehicles and attenuated by all treatments ( Figure 2A , Figure S2 , Table S7 ). Notably, microglia activation remained elevated at the CUS+Rest and CUS+US timepoints in discordance with behavioral observations while the synbiotic ubiquitously attenuated activated microglia at all time points. A similar analysis was done in the amygdala, hippocampus and nucleus of the solitary tract (NTS) of the cerebellum. In the amygdala, a similar trend was observed with respect to timepoint and treatment, while there were few inconsistent changes in both the hippocampus and the NTS ( Figures S1 , S2 , Tables S8 – 13 ). To test if microglial activation could be reflected by transcriptional control of immune genes, quantification of key sterile immune factors, canonical activators and immune cell recruitment genes in the respective brain regions was conducted. Following CUS and CUS+US, there was an upregulation of almost all immune-associated genes including toll-like receptor (Tlr)4, caspase 1 (Casp1), Il1β, Ccl2, Ccl5, intercellular adhesion molecular (Icam) and vascular adhesion molecule (Vcam) in the PFC ( Figure 2B , Table S4 ). At the gene level, the synbiotic, but not the probiotic or BDPP, robustly attenuated inflammatory gene expression profiles. Importantly, the most consistent and drastically impacted factors by both stress and the synbiotic, paralleling the behavioral response, were the immune cell recruitment genes, Ccl2, Ccl5 and Icam. In the PFC, their expression in response to stress and treatment mirrored the behavioral response with the synbiotic attenuating their stress-induced elevation. This was similarly true for the gene and cytokine expression in the PFC, hippocampus, cortex and cerebellum, however less defined compared to the PFC ( Figure 2C ). This suggests that the recruitment of peripheral immune cells into the brain may be more important than microglia activation for driving the gut-microbiota stress-induced behavioral responses.

The gut microbiome is sensitive to both stress and treatment and plays a critical role in directing the immune response along the gut-brain-axis (2). To determine how the synbiotic modulates the gut microbiome’s response to the stress protocol, 16S metagenomic sequencing was conducted. The stress protocol caused a significant shift in both the Fisher alpha (within sample variance) and beta diversities (between sample variance) in the gut microbiome of stressed vehicle controls ( Figure 3A ). The stress-induced variations in alpha and beta diversities were normalized by BDPP ( Figure 3B ) and the synbiotic ( Figure 3D ), but not probiotic treatment ( Figure 3C ). Based on these results, the chronic stress-induced variation following CUS can be interpreted as a loss of microbial diversity, which partially recovered after a period of rest. Note, that the alpha diversity was unchanged in unstressed treated controls; however all treatments elicited a similar shift in the beta diversity ( Figure S3A ). Likewise, treatment affected beta diversity following US and CUS ( Figures S3B, C ), while alpha diversity varied only following CUS due to an increase in diversity observed with probiotic and synbiotic treatment ( Figure S3C ). Interestingly, treatment had no effect on the alpha or beta diversities following CUS+Rest or CUS+US ( Figures S3D, E ). Taken together, the synbiotic was the most effective at normalizing the changes in gut microbiota diversity due to chronic and recurrent stress.

Quantification of the individual gut microbiome populations confirmed the consistent beneficial role of the synbiotic. There were no variations in either the phyla ( Figure S3F ) or genera ( Figure 3E ) of non-stressed treated controls due to either stress or treatment; however, both treatment and time elicited variations in the genera of stressed mice ( Figures 3F–I , S5 ). The most drastic variation following CUS in the stressed vehicles was an increase in Dubosiella spp., which recovered after the rest period but failed to increase upon CUS+US ( Figure S4 ). Following CUS, both probiotic and synbiotic treatments effectively reduced Dubosiella spp. levels ( Figure S5 ), while the synbiotic attenuated Dubosiella levels at all time points ( Figure S4 ). Faecalibacterium spp. is related to Dubosiella spp. yet there were no variations in Faecalibacterium spp. in stressed vehicle controls at US or CUS in contrast to the increase observed following CUS+Rest and CUS+US ( Figure S4 ). This trend was paralleled by all treatment groups ( Figure S4 ); however, the synbiotic elicited a proportionally larger increase in Faecalibacterium spp. than the other treatments ( Figures 3F , S5 ). Akkermansia spp. is another important genera in terms of metabolite production, but there were no variations due to stress in vehicle controls at any timepoint ( Figure S4 ). Importantly, the synbiotic alone elicited a positive increase in Akkermansia spp. following CUS compared to stressed vehicle controls ( Figures 3F , S5 ). Finally, the probiotic and synbiotic elicited a stark increase (and BDPP a decrease) in Lactococcus spp. at all time points, compared to the non-stressed vehicle controls ( Figure S4 ), yet there were no variations in Lactococcus spp. due to stress or timepoint in vehicle controls. Overall, the probiotic and synbiotic treatments normalized the major genera variations induced by stress in vehicle controls while the synbiotic elicited some unique changes in genera populations that play important roles in managing barrier integrity.

To understand how variations in the gut microbiota, and their interactions, may relate to the behavioral response in the context of chronic and recurrent stress, a MARS model was performed with behavior as the response and all available microbiome genera and conditions (i.e. time, stress or treatment) as predictors. In this context, the MARS algorithm takes into account the variations of all the gut microbiota populations, group and timepoint effect and the combined behavioral data to unbiasedly predict which groups of bacteria influencing the overall phenotype of chronic-stress induced behavioral impairment. The use of the MARS approach allows for the systematic identification of primary and secondary levels of interactions between gut microbiome populations, which optimize the prediction of stress-induced behavioral responses. It also controls the individual covariates’ influence (i.e. time, stress or treatment) on the behavioral response based on the microbiome data. Depressive-like behavior associated with variations in Lactobacillus spp., Ruminoclostridium spp., and secondary interactions of Lactobacillus spp. with Faecalibaculum, Blautia, or Bifidobacterium spp. ( Figure 3J ). While the log10 transformed MARS coefficients appear very small, these interactions remain significant considering the extensive size of the model with multiple covariate influences. With a larger coefficient, the interaction of stress with either probiotic or synbiotic imposed a positive effect on the depressive-like behavior ( Figure 3J ). In contrast, anxiety-like behavior did not associate to any specific genera ( Figure 3K ), but synbiotic and probiotic treatment did associate to anxiety-like behavior. Additionally, the interaction of the synbiotic with stress or the CUS timepoint both associated to lower anxiety demonstrating the beneficial effect of the synbiotic. These results indicate that probiotics and synbiotics, but not BDPP, alter the gut microbiota in a manner correlating to attenuated stress-induced anxiety- and depressive-like behavior, while depressive-like behavior is more sensitive to specific gut microbiota genera including Lactobacillus spp. Further studies will be required to understand the direct relationship between genera specially altered by the synbiotic and the direct psychological consequences.

Chronic stress and variations in the gut microbiota are intimately related to barrier immunity and programming of peripheral cellular immune responses. Further, the behavioral responses to stress and synbiotic treatment were more sensitive to the recruitment of peripheral immune cells to the brain than neuroinflammation. To characterize the peripheral innate and adaptive immune cell variations in response to stress and treatment, immunophenotyping was conducted with CyTOF. In the ileum, viSNE analysis gated by the major cell populations ( Figure 4A ) and the corresponding quantification of the major cell types ( Figure 4B ) revealed that neither treatment nor time elicited significant variations in the major ileal immune cell populations. Note that the viSNE plots were calculated for each timepoint so direct comparisons of the phenotypic islands cannot be compared between timepoints. Continuing broad characterization of the ileal immunological response, the gene and cytokine expression of key inflammatory components were determined. Following CUS and CUS+US, but not US or CUS+Rest, there was a strong upregulation of kynurenine and the proinflammatory cytokines in the ileum including IL-17A, IL-1β and IL-6, an effect rescued most strongly by the synbiotic ( Figure 4C , Table S5 ). There was a correlating downregulation of the anti-inflammatory cytokine IL-10 in control mice in response to stress, while the synbiotic ubiquitously up-regulated IL-10 release. Notably, unlike the major cell populations, cytokine expression mirrored the stress-induced behavioral responses. A similar trend was observed for gene expression. US invoked few variations in immunological gene expression, but CUS and CUS+US strongly upregulated Rorγt and Cyp1a1 in stressed vehicle controls, which were both downregulated exclusively by the synbiotic following CUS, but not CUS+US ( Figure 4D , Table S6 ). The gene expression ratio of Rorγt to Foxp3, reflecting broadly the Th17 to Treg ratio, was upregulated by CUS and CUS+US in the stressed vehicle controls, and rescued by the synbiotic following CUS and all treatments following CUS+US ( Figure S6A ). This trend was similar to the ratio of IL-17A to IL-10 cytokine release ( Figure S6B ) indicating that the variations in proinflammatory gene expression elicited by stress and beneficially affected by the synbiotic may be due to the regulation of Th17 and Treg cell differentiation.

Enumeration of Tregs and Th17 cells and their activated states was determined using the CyTOF data. In stressed vehicle controls, ileal Tregs were reduced following CUS and CUS+US, an effect ubiquitously upregulated by the synbiotic, CUS+US by BDPP and CUS+Rest and CUS+US by probiotics ( Figure 4E ). Interestingly, CTLA4⁺ Tregs (i.e. activated Tregs) were correspondingly decreased following CUS and CUS+US in stressed vehicle controls and elevated at all timepoints exclusively by the synbiotic ( Figure 4E ). A different trend was observed for the Th17 cells. Th17 cell numbers were elevated following CUS and CUS+US in the stressed vehicle controls, but only downregulated by the synbiotic in the CUS+US group ( Figure 4F ). Likewise, there were no significant changes in the IL-17A⁺ICOS⁺ Th17 cells (i.e. activated Th17) with respect to time or treatment. Concatenating this, the ratio of inactivated Th17/Treg cells in the ileum was significantly elevated in stressed vehicle controls following CUS and CUS+US, compared to unstressed vehicle controls. This stress-induced elevation was ubiquitously rescued only by the synbiotic at all time points (Fig S6c). A similar trend was observed for the ratio of “activated” Th17/Treg, except the synbiotic exclusively reduced the activated Th17/Treg ratio except at the CUS+US timepoint where probiotics also elicited a beneficial effect ( Figure S6C ). These trends were unique to the Treg and Th17 cells, as the other effector lymphocytes in the ileum including CD8 T cells (IFNγ+ and TNFα+), Th1 and Th2 cells showed inconsistent variations with respect to time and treatment ( Figure S7 , panels within).

To understand which innate immune cells may be activating the lymphocytes in the context of gut microbiota-inflammatory-stress interactions, several APCs were also assessed in the ileum. The quantities of monocytes, macrophages, dendritic cells (DC1, DC2, pDC), natural killer (NK) cells and neutrophils showed little variation with respect to either time or treatment ( Figure S7 , panels within). As reported in other studies, there was an increase in MHCII⁺ monocytes in the ileum of stressed vehicle controls following CUS and CUS+Rest, which was reduced by all treatments ( Figure S7C ); however, there was no effect following CUS+US. Interestingly, ILC, ILC1, ILC2 ( Figure S7 ) and ILC3 ( Figure 4G ) populations elicited few variations over time and treatment; however, the IL-22⁺ (NCR⁺) and IL-17A⁺ (NCR^-) producing ILC3 cells showed trends paralleling the behavioral responses. IL-22 producing ILC3s were significantly downregulated in response to stress in vehicle controls, while synbiotic treatment exclusively upregulated their expression ( Figure S7J ). Likewise, synbiotic treatment downregulated the stress-induced increase in IL-17A producing ILC3 cells ( Figure S7K ). Overall, the ratio of IL-22/IL-17A producing ILC3s was reduced in response to stress and rescued exclusively by the synbiotic at both the CUS and CUS+US timepoints ( Figure 4H ). This is a significant result as IL-22 and IL-17A producing ILC3s are, in part, responsible for the differentiation of either Treg or Th17, respectively (45). It also suggests that manipulation of the gut microbiota with a synbiotic may reprogram the barrier immune architecture towards an NCR+ ILC3-Treg dominating phenotype that consequently prevents the stress-induced inflammatory milieu associated with chronic and recurrent stress.

To expand on the immunophenotyping of the ileum and gather a full understanding of the microbiota-peripheral-brain immune crosstalk, innate and adaptive immune cells in both the spleen ( Figure 5A , Figure S8A , Tables S14 , 15 ) and PBMCs ( Figure 5B , Figure S8D , Tables S16 , 17 ) were also immunophenotyped. In the spleen, stress-induced alterations in the Th17/Treg activated cell ratio was similar to the ileum, except the synbiotic elicited a less drastic impact ( Figure 5C ), which can be attributed to variations in the Treg population ( Figure S9G ). Variations in the Th17/Treg ratio were not observed among the PBMCs ( Figure 5F ). There was also an increase in MHCII⁺ monocytes in the spleen following CUS in the stressed vehicle controls ( Figure 5D ) and a trending decrease of migratory CD103⁺ DC1 cells, where the latter was rescued by the synbiotic following CUS and CUS+Rest ( Figure 5E ). There was a similar increase in activated monocytes among PBMCs following CUS in the stressed vehicle control, an effect rescued by synbiotics alone ( Figure 5G ); but no alterations due to stress or treatment were observed in the migratory DC1 population ( Figure 5H ).

Gene and cytokine expression previously shown to be altered by stress and synbiotic treatment in the periphery (28), showed less consistent variations in the context of chronic and recurrent stress. There was a stark increase in Rorγt, Cyp1a1 and Tlr4 in the splenic stressed vehicle controls following CUS, all rescued by synbiotic treatment ( Figure S8B ). Importantly, these variations were not observed in the CUS+US group uncoupling the inflammatory response in the spleen from the behavioral response. A similar trend was observed for splenic cytokine expression ( Figure S8C ). In the serum, the proinflammatory cytokines IL-1β and IL-17A were upregulated following both CUS and CUS+US in the stressed vehicle controls, and downregulated by the synbiotic ( Figure S8E ). Notably, IL-17A expression remained elevated following rest (CUS+Rest) in the stressed vehicle controls, also uncoupling it from the behavioral phenotype. Importantly, unlike in the ileum, the variations in immune cells observed in the PBMCs and spleen generally did not associate with the chronic and recurrent stress-induced behavioral impairment, yet a trend towards a chronic-stress induced peripheral immune cell phenotype may be apparent. As such, in both the spleen ( Figure S9 ) and PBMCs ( Figure S10 ), variations in minor cell populations were observed, yet not obviously correlated to the behavioral phenotypes.

To systematically dissect which minor cell populations were varying with respect to treatment and/or time, the MARS model was repeated, using an additive model, with the cell frequencies as the predictors to the timepoint, treatment and stress as responses. As these data were collected from a separate cohort of animals not used for the above behavioral analyses, exploration of a link between cell frequencies and behavior was not attempted. This analysis will determine which of the cell frequencies are associated to the study conditions, i.e. day, stress or treatment. In the ileum, the Th17/Treg ratio was positively associated to the control and BDPP treatments, while negatively associated to probiotic and synbiotic, defining the beneficial effect specifically of the probiotic and synbiotic on this lymphocyte ratio within the context of chronic and recurrent stress ( Figure 5I ). A similar trend was observed for the activated CTLA4+ activated Tregs. Regarding the variation associated to timepoint (CUS, CUS+Rest or CUS+US), the monocytes, macrophages and dendritic cells all displayed significant associations. This apparent dichotomy between cell types responsive to timepoint (mostly APCs), cell types associated to treatment (lymphocytes and ILCs) and cell types varying with stress demonstrate a clear separation of regulatory mechanisms to be considered in the ileum of the current stress and treatment model.

A similar separation was observed in the spleen and PBMCS. In the spleen, there was a strong association of the immune cell subsets with the timepoint, while only macrophages weakly associated with treatment ( Figure 5J ). In addition, and in support of the cell frequency quantification, MHCII⁺ monocytes were highly associated with the overall stress condition, whereas MHCII⁺ macrophages more strongly associated to the individual timepoints. Finally, cell subsets in stressed animals significantly varied with respect to timepoint in the PBMCs, not treatment and only weakly with stress ( Figure 5K ). The strongest association in PBMCs was the activated MHCII⁺ dendritic cells and the Tregs, which were both responsive to timepoints associated with stress. This analysis revealed that there is an immune cell tissue-dependency in the context of chronic and recurrent stress, with the ileum having a considerable association to both time and treatment and in the spleen and PBMCs, the associations are primarily with the timepoint.

The peripheral, especially ileal, immune response is an important predictor of chronic and recurrent-stress induced behavioral impairment; however, neuroinflammatory markers in limbic brain regions are common signatures of stress-induced depressive- and anxiety-like behaviors. To dissect which peripheral and/or brain-derived immune factors are associated to the chronic and recurrent stress-induced behavioral responses, linear and multivariate interactions between immune factors and the behavioral response were made. Data was only included which was collected from the same animal, allowing only the 16S sequencing, behavioral, RNA and protein enumerations to be included. Linear associations between the gene and cytokine factors with the behavior response were conducted, with the gene and cytokines acting as the response variable and the behavior as the main predictor. In this model, the stress, timepoint and treatment were all treated as covariates giving an overall indication of the how the behavioral response associates with single tissue-specific immune genes and cytokines. Note, that the effect size represents the strength of the association while correlation represents the overall fit of the model in the context of all the covariates. In both cases where the forced swim test ( Figure 6A ) or the open field test ( Figure 6B ) were used as the predictor, one striking observation was that there were more peripheral immune associations compared to brain-derived ones. Predicting depressive-like behavior, only Il1β gene expression in the PFC and cerebellar Ccl5 elicited a positive effect size, with the PFC response having a higher correlation than the cerebellum, indicating that the entire model fit well for the behavioral prediction. Splenic IL-6 (gene and cytokine) and Cyp1a1 also displayed positive effect sizes with splenic Cyp1a1 having a strong correlation demonstrating the importance of AHR pathway signaling. The remaining seven associations occurred in the ileum with the proinflammatory markers kynurenine, IL-17, IL-1β and IL-6 having positive and Foxp3 and Gata3 gene expression negative effect sizes ( Figure 6A ). Predicting anxiety-like behavior, which has an opposite direction of beneficial effect compared to the forced swim test (i.e. “lower” value in the open field test means less anxiety-like behavior), cerebellar IL-6 had a positive effect size and good correlation while Icam, Il10 and Il1β all had negative effect sizes and weak correlations ( Figure 6B ). Similar to the forced swim test data, splenic Cyp1a1 had a negative effect size and Gata3 gene expression a positive effect size and strong correlation to the model. Also similar to the depressive-like behavior, was the dominant representation of ileal immune associations. Again, Foxp3 had a strong response to anxiety-like behavior in the ileum with an opposing effect size compared to the other proinflammatory markers such as kynurenine, Nfkb, Tlr2, Il1β and Il10 ( Figure 6B ). Overall, the linear model confirms that the peripheral, namely ileal, immune responses are better associated than the limbic brain-derived immune factors to the chronic and recurrent stress-induced behavioral response within the context of gut microbiota modifying prebiotics, probiotics and synbiotic suggesting that a direct relationship may exist.

Linear associations give valuable insights to the impact of the behavioral response on a single gene or cytokine predictor; however, this model is limited as the genes and cytokines are assumed to be independent, which could contribute to the positive association or explain the low associative coefficient. To detangle these covariates and predict which individual factors across tissue, treatment and time may be synergizing in response to the behavioral output, the MARS approach was implemented. Similar to the linear model, behavior [either depressive-like ( Figure 6C ) or anxiety-like ( Figure 6D )] was used as a response variable, while all available genes, cytokines and sample characteristics (timepoint, stress, and treatment) were used as predictors. Using the generalized R-squared as criterion, allowing a maximum of three interactions in a single term provided a better model fit than when a maximum of two interactions were allowed (see details in materials and methods section). Using this analysis, primary, secondary and tertiary interactions can be identified. Similar to the linear associations, the ileum had the highest predictive power, with genes (Rorγt, IL-23R and Nfκb) and cytokines (IL-17A and IL-1β) having significant primary associations in both behaviors tested. For the spleen, gene expression of Il6, Il1β and Nfκb likewise associated with behavior. In the brain, IL-1β in the PFC associated with depressive-like behavior whereas Vcam, Tlr2 and IL-1β in the cortex and Il-10 in the cerebellum also had significant primary effects.

The tertiary interactions bring a novel understanding to the capacity of the entire gut-brain-axis to drive chronic and recurrent stress-induced depressive-like behaviors synergizing across the tissues investigated. Interestingly, the triple association between ileal Foxp3, Rorγt and hippocampal Ccl5 strongly associated with depressive-like behavior ( Figure 6C ) linking the Th17/Treg ratio to the recruitment of peripheral immune cells to the hippocampus. Likewise, ileal Rorγt, cerebellar Ccl5 and hippocampal Ccl5 were a strong predictor of depressive-like behavior as was ileal Cyp1a1, splenic Il10 and PFC IL-1β. Ileal Nfkb expression also coordinated with hippocampal Ccl5 and IL-6 production to predict anxiety-like behavior ( Figure 6D ). A splenic increase in Il1β and Tlr4 with cortical Vcam expression strongly associated with the anxiety-like behavior while splenic IL-10 associated with IL-1β from both the PFC and cortex, predicted depressive-like behavior. These and the other associations illustrate how the gut-brain-axis may functionally coordinate with the limbic brain regions to predict stress-induced depressive- and anxiety-like behavior outlining the importance of the gut microbiota and barrier immunity in driving these responses. These correlations set the foundation for further studies to directly validate the impact of these singular and synergistic associations using transgenic mice and specifically immuno-compromised mice.