Social isolation alters behavior, the gut-immune-brain axis, and neurochemical circuits in male and female prairie voles.

The absence of social support, or social isolation, can be stressful, leading to a suite of physical and psychological health issues. Growing evidence suggests that disruption of the gut-immune-brain axis plays a crucial role in the negative outcomes seen from social isolation stress. However, the mechanisms remain largely unknown. The socially monogamous prairie vole (Microtus ochrogaster) has been validated as a useful model for studying negative effects of social isolation on the brain and behaviors, yet how the gut microbiome and central immune system are altered in isolated prairie voles are still unknown. Here, we utilized this social rodent to examine how social isolation stress alters the gut-immune-brain axis and relevant behaviors. Adult male and female prairie voles (n = 48 per sex) experienced social isolation or were cohoused with a same-sex cagemate (control) for six weeks. Thereafter, their social and anxiety-like behaviors, neuronal circuit activation, neurochemical expression, and microgliosis in key brain regions, as well as gut microbiome alterations from the isolation treatment were examined. Social isolation increased anxiety-like behaviors and impaired social affiliation. Isolation also resulted in sex- and brain region-specific alterations in neuronal activation, neurochemical expression, and microgliosis. Further, social isolation resulted in alterations to the gut microbiome that were correlated with key brain and behavioral measures. Our data suggest that social isolation alters the gut-immune-brain axis in a sex-dependent manner and that gut microbes, central glial cells, and neurochemical systems may play a critical, integrative role in mediating negative outcomes from social isolation.

Introduction

Humans, who have evolved to be highly social, rely on communities for support throughout life. In contrast, a lack of social bonds and support from others can be particularly detrimental to one's health. Social isolation, or the lack of social interactions, is becoming an exponentially widespread issue with a wide variety of negative outcomes (Nicholson, 2012). For example, social isolation and loneliness results in increased rates of depression and anxiety in humans (Domènech-Abella et al., 2019; Ge et al., 2017; Wang et al., 2018). Perceived social isolation can also greatly influence one's inflammation levels, suggesting a strong tie between social bonds, the immune system, and physical health outcomes (Hawkley and Capitanio, 2015; Yang et al., 2013). In fact, perceived isolation is considered a strong predictor of mortality, in line with other commonly discussed clinical issues, such as obesity (Holt-Lunstad et al., 2015; Laugesen et al., 2018; Pantell et al., 2013). As social isolation can both activate the hypothalamus-pituitary-adrenal (HPA) axis (Hawkley et al., 2012) and alter immune function in the body (Eisenberger et al., 2017; Gądek-Michalska et al., 2017), more detailed investigations of how social isolation stress, the brain, and the immune system are interrelated are needed. Further, as the majority of the immune system resides in the gut (Vighi et al., 2008), a focus on how gut microbes may also be intricately tied to stressor-induced alterations is becoming increasingly prevalent (Foster et al., 2017). These microbes within the gut can alter brain function through a variety of pathways, including vagus nerve activation, metabolite circulation, and immune alterations; collectively, these bidirectional pathways are referred to as the microbiota-gut-brain axis (Cryan and O'Mahony, 2011; Mayer et al., 2015). Given the dramatic rise in social stressors in today's societies (Yang et al., 2013), a better understanding of how social isolation impacts the gut-immune-brain axis and results in atypical behaviors is needed.

Social isolation-induced effects on the brain and behaviors have been well documented in animal models. Social isolation and resulting anxiety have been shown to alter a wide variety of neurochemical systems in the brain, including oxytocin (OT), serotonin (5-HT), and corticotropin releasing factor (CRF) systems in various brain regions known to be involved in social and stress-related behaviors (Pan et al., 2009; Pournajafi-Nazarloo et al., 2011). This isolation-induced anxiety can also be detrimental to the immune system (Cohen et al., 1997; Glaser et al., 1985). For example, in rats, increased anxiety-like behavior by social isolation is linked to immunosuppression (Cruces et al., 2014), and social isolation in early life increases depressive-like behaviors via microglial activation in the hippocampus, which can be reversed if the microglial activation is blocked (Wang et al., 2017). As microglia are deemed the resident immune cells of the central nervous system (Lenz and Nelson, 2018), these data indicate both central and peripheral immune alterations by social isolation. Further, social isolation can have a differential impact on males and females, and such sex differences can also be tied to differences in inflammatory activation (Hermes et al., 2006; Weiss et al., 2004). A few experiments have started examining the impacts of social stressors on the gut microbiome using animal models. For example, a recent study found that adolescent social isolation significantly alters gut microbiome composition, which correlates with an increase in depressive-like behavior in rats (Dunphy-Doherty et al., 2018). Other social stressors, such as social defeat, can alter the gut microbiome and neuroinflammation in mice as well (McKim et al., 2016; Yang et al., 2017). The large body of research demonstrating the profoundly negative effects of social isolation on neurochemical systems in the brain, neuroinflammation, and behaviors in traditional rodent models suggests a critical link among these systems.

The socially monogamous prairie vole (Microtus ochrogaster) is a rodent species that naturally displays a wide variety of social behaviors, and thus has been established as an alternate, translational animal model for studying social behavior (Young et al., 2011). As prairie voles are highly social in nature, the absence of social bonds can be particularly detrimental. Indeed, social isolation has been shown to increase anxiety- and depressive-like behaviors and to alter neurochemical systems in the brain in prairie voles (Grippo et al, 2007b, 2008, 2009, 2015; Lieberwirth et al., 2012; Scotti et al., 2015). Interestingly, social isolation can also alter the peripheral immune response to a challenge in prairie voles (Scotti et al., 2015). Yet, no studies have examined how social isolation alters central immune function or gut microbiota composition in prairie voles. The prairie vole gut has been previously shown to consist of common predominant taxa at the phylum level, such as Firmicutes and Bacteroidetes (Curtis et al., 2018). The prairie vole microbiome also contains common microbes deemed beneficial for human health, such as Lactobacillus species (Assefa et al., 2015). We also used metagenomic sequencing to better understand the prairie vole microbiome and found it to contain novel, dominant strains, suggesting that the prairie vole gut composition is quite unique from rats and mice (Donovan et al., 2020). As the microbiome is highly intertwined with the immune system (Galley and Bailey, 2014; Geuking et al., 2014; Palm et al., 2015), social stress can alter the gut microbiome in traditional rodent models (Bailey et al., 2011; Bharwani et al., 2016; Burokas et al., 2017; Partrick et al., 2018), and specific bacterial strains can alter various inflammatory markers in the body (Mackos et al., 2013; Thomas et al., 2012), a closer examination of how microglia and gut microbiota are altered from social isolation using the prairie vole model may reveal a link between the gut-immune-axis and behaviors.

In this study, we focus on microgliosis as our primary immune measurement, as it has been heavily implicated in inflammatory responses to social stressors (Krügel et al., 2014; McKim et al., 2016; Wang et al., 2017). There have been exciting developments demonstrating the link between microbiota alterations and effects on microglia activity and resulting neuroinflammation (D'Mello and Swain, 2017; Thion et al., 2018). In fact, microbiota are necessary for proper microglia maturation, functioning, and development (Erny et al., 2015). Furthermore, despite well-known sex differences in animals' bodily functions and in their responses to altered environments (including social environments), the majority of previous studies only focus on one sex of experimental animals. Here, we used our established six-week social isolation paradigm (Lieberwirth et al., 2012; Pan et al., 2009) to systematically examine effects on the gut microbiome, neurochemical and neuroinflammatory alterations in the brain, and behaviors in male and female prairie voles. This study is the first to examine both microglial alterations and the gut microbiome post-isolation with behavioral, neuronal, and neurochemical measurements. An assessment of how microgliosis in key brain regions and the gut microbiome are altered from a social stressor using a translational animal model will enhance our knowledge of their involvement in social behaviors (Sherwin et al., 2019). We hypothesized that social isolation would increase anxiety-like and alter social behaviors in prairie voles, and these behavioral changes would be associated with altered neurochemical circuit activation, increased microgliosis in key brain regions for social behavior, and differential microbiome composition. In terms of microbiome diversity, we hypothesized that there would be no significant differences between the two groups at baseline, but there would be significant changes to microbiome diversity across cohoused and isolated subjects post-treatment.

Methods & materials

Subjects and social isolation paradigm

Subjects were male and female prairie voles (M. ochrogaster) captive-bred at Florida State University. Subjects were weaned on postnatal day 21 and housed in Plexiglas cages (20 × 25 × 45 cm) with a same-sex conspecific. All cages contained cedar chip bedding with food and water provided ad libitum. All subjects were kept at 20 °C under a 14:10 h light:dark cycle (lights on at 0700). At the time of isolation, subjects had reached adulthood (>90 days of age) and were sexually naïve. The isolation procedure has been established in our previous study (Lieberwirth et al., 2012; Pan et al., 2009). Briefly, sexually naïve prairie voles were randomly placed into a clean cage containing cedar chip bedding with food and water ad libitum either alone or with their same sex cagemate. The isolation procedure lasted for 6 wks, and behavioral testing was performed immediately after. All procedures were approved by the Institutional Animal Care and Use Committee at Florida State University and were in accordance with the guidelines set forth by the National Institutes of Health.

Behavioral testing

The elevated plus maze (EPM) test has been established and validated in our previous vole studies to examine anxiety-like behaviors (Smith and Wang, 2014). The apparatus is elevated 45 cm off the ground and consists of two open (35 × 6.5 cm) and two closed arms (35 × 5 × 15 (H) cm) that cross in the middle. After 6 wks of housing treatment, subjects were placed onto the center of the maze, and subject's behaviors were recorded for 5 min using Active Webcam software. Behaviors quantified, including the duration and frequency in the open arms, closed arms, and in the center of the maze, were scored by an observer blind to treatment via J-Watcher. The percentage of time in the open arms and locomotor activity (total entries) were also calculated.

The social affiliation (SA) test has also been established in our previous vole study to examine social affiliation behavior (Pan et al., 2009). The apparatus consists of two polycarbonate chambers (13 × 18 × 29 (H) cm) connected by a hollow tube (7.5 × 16 cm). A same-sex, unfamiliar stimulus animal at a similar age and size as the subject was loosely tethered in one chamber (counterbalanced across subjects), and the subject was placed into the empty chamber to freely roam the apparatus. Each cage contained fresh cedar chip bedding. The SA test lasted for 30 min during which subjects’ behaviors were recorded using Active Webcam software and were subsequently quantified by an observer blind to treatment via J-Watcher. Behaviors quantified included duration and frequency in the conspecific cage, empty cage, and connecting tube. Additional behavioral quantifications included duration and frequency of specific behaviors during the test (Table S1).

Brain tissue preparation

All subjects that underwent EPM testing were rapidly decapitated 1 h after the start of behavioral testing. Brains were immediately extracted, placed on dry ice, and were then stored at −80 °C until processing. Tissue punches were collected from coronal sections (200 μm thick) from the prefrontal cortex (PFC), amygdala (AMY), nucleus accumbens (NAcc), paraventricular nucleus of the hypothalamus (PVN), and dentate gyrus (DG) (four sections per region per subject). Tissue punches were stored at −80 °C until subsequent protein extraction. All subjects that underwent SA testing were perfused 1 h after the start of behavioral testing. Subjects were anesthetized and transcardially perfused via 0.9% saline followed by 4% paraformaldehyde solution. Brains were then collected and fixed in 4% paraformaldehyde for 2 h until placed in 30% sucrose in phosphate buffer (PB) at 4 °C for storage. Tissue slices were cut into 40 μm coronal sections via sliding microtome. Each subject only underwent one behavioral test (see Experimental Design below).

Protein extraction

Protein extractions were performed using established methods (Tabbaa et al., 2017a). Briefly, protein was extracted from tissue punches using Tri-Reagent according to the manufacturer instructions (Molecular Research Center, Cincinnati, OH). Protein was stored at −80 °C until western blotting.

Western blotting

Western blotting procedures were performed using our established methods (Tabbaa et al., 2017a). Briefly, 15 μg of protein were loaded into 15% sodium dodecyl sulfur (SDS) polyacrylamide gels (Bio-Rad) for electrophoresis. Proteins were run on gels at 75 V (V) for approximately 30 min and thereafter at 200 V for 1 h 20 min under refrigeration. Proteins were then transferred to nitrocellulose membranes and subsequently blocked in 5% milk or Superblock (Bio-Rad). Membranes were then incubated for 1–2 days with one of the following primary antibodies: rabbit anti-glyceraldehyde 3-phosphate dehydrogenase (GAPDH) (1:1 K, Santa Cruz), goat anti-oxytocin receptor (OTR) (1:1 K, Santa Cruz), rabbit anti-CRF (1:500, ProteinTech), rabbit anti-corticotropin releasing factor receptor 1 (CRFR1) (1:500, Novus Biologicals), rabbit anti-corticotropin releasing factor receptor 2 (CRFR2) (1:500, Novus Biologicals), brain derived neurotropic factor (BDNF) (1:500, Santa Cruz), and receptor tropomyosin receptor kinase B (TrkB) (1:1 K, Santa Cruz). Thereafter, membranes were washed in TPBS and incubated for 2 h in respective horseradish peroxidase (HRP) conjugated secondary antibodies (1:10 K, Santa Cruz). Membranes were then washed again in TPBS for 1 h. All membranes were then placed in chemiluminescence HRP substrate (SuperSignal West Dura Extended Duration Substrate, Thermo Fisher Scientific) for 10 min. All bands were visualized on the ChemiDoc MP System (Bio-Rad). Western band quantification was done through ImageLab software. All markers were normalized using GAPDH as a loading control. Data are expressed as a ratio of marker over GAPDH signal. All antibodies have been validated and used in our previous study in prairie voles (Tabbaa et al., 2017a).

Immunocytochemistry

Ionized calcium binding adaptor molecule 1 (Iba-1) and Early Growth Response protein 1 (Egr-1) staining were performed on sets of coronal sections with 200 μm intervals using our previously established method (Liu et al., 2019). Briefly, sections were rinsed in 0.1 M PBS 4 times for a total 20 min (each rinsing was kept at the rate of 4 times for a total of 20 min below) and then treated with 1% NaBH4 in 0.1 M PBS for 10 min. After rinsing in 0.1 M PBS, sections were incubated in 0.3% H2O2 in 0.1 M PBS for 20 min. After rinsing in 0.1 M PBS, sections were then incubated in 10% normal goat serum (NGS, Sigma-Aldrich, St. Luis, MO) in 0.3% Triton X-100 in 0.1 M PBS (TPBS) for 1 h at room temperature. Sections were then incubated in rabbit anti-Egr-1 (1:3 K, Cell Signaling Technology, Inc., Danvers, MA) or rabbit anti-Iba-1 (1:10 K, Wako Chemicals) in 0.3% TPBS with 2% NGS at 4 °C for 2 nights. Sections were then placed at room temperature for 1 h and then rinsed in 0.3% TPBS. Sections were then incubated in biotinylated goat antirabbit IgG (1:300, Vector Laboratories, Burlingame, CA) for 2 h at room temperature. Thereafter, sections were rinsed in 0.3% TPBS followed by 0.1 M PBS and were then incubated in the ABC Elite HRP Kit (Vector Laboratories) in 0.1 M PBS for 90 min at room temperature. After 0.1 M PBS rinses, Egr-1 and Iba-1 immunostaining was revealed using 3′-diaminobenzidine (DAB, Sigma-Aldrich). Double-label staining for Egr-1/OT, Egr-1/AVP (vasopressin), Egr-1/CRF, and Egr-1/5-HT were also performed on sets of coronal sections with 200 μm intervals using previously established methods (Liu et al., 2019). For double labeling, sections were rinsed after DAB staining in 0.1 M PBS and then incubated in 5% NGS or 5% NRS in 0.3% TPBS for 30 min. Sections were then incubated in one of the following primary antibodies at 4 °C for 2 nights: OT (ImmunoStar, 1:50 K), AVP (Millipore, 1:8 K), CRF (Peninsula Lab, 1:3 K), or 5-HT (ImmunoStar, 1:60 K). Sections then underwent the same second day protocol from the Egr-1 protocol, but were stained using SG (Vectorlabs, SK-4700) instead of DAB. After staining, all sections were mounted onto slides, dehydrated in ethanol, clarified in xylene, and cover slipped with permount. The specificity of the Iba-1 antibody has been previously validated in prairie voles (Rebuli et al., 2016) and was also validated in our lab by omission of the primary antibody, which did not lead to any specific staining.

Blood preparation and corticosterone radioimmunoassay

Blood samples (~300 μl) were collected in microcentrifuge vials containing 20 μl EDTA and immediately placed on ice. Samples were centrifuged at 6000 rpm for 15 min at 4 °C, and plasma was transferred into new tubes. Plasma was then re-centrifuged at 6000 rpm for 10 min at 4 °C. Plasma samples (1:600) were analyzed for corticosterone (CORT) via radioimmunoassay (RIA) using a commercially available kit according to manufacturer instructions (Diagnostic Products Corp., Los Angeles, CA). The kit has been previously validated in prairie vole studies (Smith et al., 2013; Smith and Wang, 2014).

Stool sample collection and DNA extraction

Stool samples were collected by placing subjects in a clean Plexiglass cage (20 × 25 × 45 cm) containing an elevated wired mesh surface (http://www.ancare.com). A sterile pad was placed underneath the wired mesh to allow for sterile stool collection. Subjects were placed on the platform and allowed to freely roam the cage for 1 h so that an adequate amount of stool samples could fall onto the sterile pad. Subjects were then removed from the apparatus, and stool samples were collected into sterile tubes. All samples were stored at −80 °C until further processing. Pre-treatment collections occurred 1 day pre-isolation treatment, and post-treatment collections occurred 1 day before behavioral testing.

DNA was prepared from frozen stool samples using the MoBio/QiaAmp PowerFecal DNA kit according to manufacturer's instructions (Qiagen-USA, Germantown, MD). Isolated DNA was quantified by absorbance at 260 nm on a Nanodrop spectrophotometer and by fluorescence using the dsDNA HS DNA Assay on a Qubit fluorometer (both instruments by ThermoFisher, Waltham, MA, USA).

16 S rDNA library construction and sequence analyses

16 S rDNA libraries were prepared according to Illumina's 16 S Metagenomic Sequencing Library Preparation reference guide (part# 15044223 B; https://support.illumina.com/downloads/16s_metagenomic_sequencing_library_preparation.html). V3–V4 region amplicon primers S-D-Bact-0341-b-S-17, 5′-CCTACGGGNGGCWGCAG-3′ and S-D-Bact-0785-a-A-21, 5′-GACTACHVGGGTATCTAATCC-3’ (Klindworth et al., 2013) were modified to include a heterogeneity spacer containing 0 to 3 random nucleotides between the 16 S sequence and the adaptor sequence, in order to create more sequence diversity and eliminate the need for phiX spike-in (Fadrosh et al., 2014). Libraries were pooled and sequenced as 2 × 300 bp reads using Illumina MiSeq v3 600 cycle reagents according to manufacturer's instructions. Sequence data is deposited at NCBI as BioProject PRJNA629975. Accession numbers for individual animal samples and raw reads can be found in Supplemental Material.

Raw paired-end 16 S rDNA sequences were trimmed for adapter sequences using the Qiime2 (Bolyen et al., 2018; Caporaso et al., 2010) cutadapt plugin (Martin, 2011) and then denoised and clustered into ASVs using the Dada2 wrapper (Callahan et al., 2016), with sequences truncated with the wrapper to a length of 280 based on quality score summaries visualized with Multi-QC (Ewels et al., 2016). Clustered sequences were then classified in the Qiime2 environment using a naïve-Bayes classifier trained off the 16s SILVA rDNA 97% similarity database (Quast et al., 2013). Classified sequences were further summarized into feature count tables and exported to allow transposition and manipulation into treatment group subsets. Further analyses were then performed outside the Qiime2 environment.

Differential representation analysis of taxa in the 16 S rDNA dataset

The feature count table was normalized using the DESeq2 (Love et al., 2014) R package to extract the mean log ratio size factors which the counts were divided against. The normalized count table was then subset into groups to allow comparisons between groups. Subset count tables were then analyzed using the DESeq2 R package's Differential Expression wrapper (Love et al., 2014) and using Benjamini-Hochberg False Discovery Rate procedure to correct for multiple comparisons (Benjamini and Hochberg, 1995). The false positive correction used the adjusted p value cutoff of 0.1. A separate method designated ‘Vole-by-Vole’ on the DESeq2-normalized count tables was also performed using Microsoft Excel. Briefly, in this method, differences in representation of each microbial taxon were determined for each individual animal pre-to post-social environment manipulation. For both methods, a filtering criteria of an average log2 fold change of ≥0.693 or ≤ −0.693 was used to select taxa of interest. Additionally, the ‘Vole-by-Vole’ pairwise method takes into account the individual variability of animals by requiring that each taxon have a supermajority of animals in the group each meet the |0.693| log2 fold change criteria while also requiring that the average log2 fold change be |0.223| in the same direction. Supermajority cutoffs were: 14/23 for the isolated group, 16/24 for the cohoused group, 7/11 for the isolated females, and 8/12 for all other sex-separated housing groups (isolated males, cohoused males, and cohoused females).

Mean relative abundance of taxa in the different sample groups was calculated using the Phyloseq R package (McMurdie and Holmes, 2013). First, samples were merged into the treatment groups of interest using phyloseq:merge_samples and counts were transformed into relative abundance counts by recursively dividing the raw counts by the total sum of counts for each treatment group through phyloseq:transform_sample_counts. Then, taxonomic hits were agglomerated into the taxa level desired using phyloseq:tax_glom. The results were plotted using the ggplot2 package after cleaning the data by pruning all taxonomic hits with a total sum of less than 0.01 (Wickham, 2011). Abundance quartiles were determined from average pre-manipulation raw taxa counts.

Microbial diversity analyses

Both the alpha and beta diversity analyses were performed using the Phyloseq R package (McMurdie and Holmes, 2013) supplemented with both the vegan (Oksanen et al., 2013) and ape (Paradis et al., 2004) packages. The Chao1 (Chao, 1984) estimate for richness, the Shannon (Lemos et al., 2011; Ludwig and Reynolds, 1988; Magurran, 2004; Shannon and Weaver, 1949) effective number of species (ENS), and the Inverse Simpson (1949) diversity metrics were used for alpha diversity comparisons. Chao1 is a non-parametric estimator of richness dependent on singletons and doubletons present in the un-normalized count data (Chao, 1984), and both the Shannon ENS and Inverse Simpson's are considered true diversity metrics based on the Shannon's (Lemos et al., 2011; Ludwig and Reynolds, 1988; Magurran, 2004; Shannon and Weaver, 1949) and Simpson's (Simpson, 1949) diversity estimations of un-normalized count data, respectively. Since the Shannon's diversity metric weighs species richness over evenness and Simpson's weighs species evenness over richness, both were used to balance the analyses. Significant differences were defined as those with a p value ≤ 0.05 in two-tailed Mann-Whitney tests for two independent samples.

Beta diversity (differences between samples) was analyzed by several methods: the Bray-Curtis dissimilarity, which quantifies the difference in taxa abundance without regard to phylogenetic relationships between taxa (Bray and Curtis, 1957); the Unweighted Unique fraction metric (UniFrac), which plots all taxa in the compared samples on a phylogenetic tree and determines shared branch lengths between them (Lozupone and Knight, 2005); and the Weighted UniFrac, which is similar to UniFrac, but gives greater weight to branch lengths corresponding to abundant lineages (Lozupone et al., 2007). To determine which variables have a significant effect on variation in beta diversity, a permutational multivariate analysis of variance (PERMANOVA) test was performed using vegan:adonis2. Adonis2 applies a linear model to a given distance matrix using a sum of squares method for calculating R2 values to determine goodness of fit, and an F-test to determine significance (Oksanen et al., 2015). The contribution of variables to total variation between samples at the species level was also analyzed using CAP ordination. Multi-dimensional ordination plots from Constrained Analysis of Principal coordinates (CAP) (Anderson and Willis, 2003) were created with the vegan:capscale R function to plot the beta-diversity and to test hypotheses on the contribution of variables of interest to total variation (Morelan et al., 2019; Oksanen et al., 2015; Xia et al., 2018). Given the large individual variation across voles in microbiome composition, we compared post-treatment to pre-treatment baseline stool samples for each subject, and subjects in the cohoused condition were analyzed individually.

Data quantification and analyses

All behavioral data were scored via JWatcher software program (V1.0, Macquarie University and UCLA: http://www.jwatcher.ucla.edu/) All behavioral, neurochemical, and immunolabeling data were analyzed via two-way ANOVA (SEX x TREATMENT). Significant sex-by-treatment interactions were further analyzed by the Student-Newman-Keuls posthoc tests. Egr-1-labeled cells were quantified bilaterally (using Zeiss Axioskop 2, Carl Zeiss Microscopy Ltd, Germany) in the entire NAcc, dorsal raphe (DR), and in a specified region of interest in the DG. Iba-1-labeled cells were quantified bilaterally in the entire PVN and in a specified region of interest in the NAcc, DG, and DR. For double-label experiments, all single-labeled Egr-1, single-labeled neurochemical markers, and double-labeled cells were counted in the entire PVN and DR. For western blots, all marker densities were determined using ImageLab software and normalized using GAPDH as a loading control. Western data were analyzed as a ratio of marker over GAPDH signal. Correlations were performed in the Phyloseq R package using the Spearman correlation method (Salkind, 2012) which treats data in an ordinal fashion. Correlations of brain region and behavioral phenotypes with microbial taxa were performed both with ranked post-treatment abundance and also with ranked pre-to post-treatment log2 fold change. Correlation significance was determined with the Spearman rho test using an alpha of 0.05 in the Phyloseq package. To visualize potential patterns of interest, heatmaps were created of the correlation data with the help of the ggplot2 (Wickham, 2011) and reshape 2 (Wickham, 2007) packages.

Experimental designs

Experiment 1: does social isolation alter social affiliation and neuronal circuit activation in the brain of male and female prairie voles?

Adult male and female prairie voles were randomly assigned into an experimental group that was either continuously housed with their cage mate (control, n = 12 per sex) or singly housed (isolation, n = 12 per sex) for 6 wks (Lieberwirth et al., 2012; Pan et al., 2009). Thereafter, they underwent the 30-min SA test to measure affiliative behaviors. Thirty mins following the SA test, trunk blood was collected and subjects were perfused. Their brain sections at 200 μm intervals were processed for immunohistochemistry. Single labeling for Egr-1, an immediate early gene and indirect marker of neuronal activation, and double labeling for Egr-1/OT, Egr-1/AVP, Egr-1/CRF, and Egr-1/5-HT were performed due to their known involvement in social and anxiety-like behaviors (Gobrogge and Wang, 2015; Young et al., 2011). Plasma samples were processed for CORT measurement.

Experiment 2: does social isolation alter anxiety-like behaviors and neurochemical expression in the brain of male and female prairie voles?

Adult male and female prairie voles from the 6-wk social isolation (n = 12 per sex) or cohousing (n = 12 per sex) treatment groups, as described for Experiment 1, underwent the 5-min EPM test for anxiety-like behaviors. Subjects were sacrificed 1 h after the EPM test, brains and trunk blood were collected, and tissue punches from the selected brain areas were processed for Western Blotting to quantify the expression of various neurochemicals and receptors including BDNF, CRF, CRF, CRFR1, CRFR2, OTR, and TrkB in the PVN, PFC, NAcc, AMY, and DG. Plasma samples were processed for CORT measurement.

Experiment 3: does social isolation alter microgliosis in key brain regions and the gut microbiome in a sex-dependent manner?

Stool samples were collected from all subjects in Experiments 1 and 2 the day prior to the 6-wk treatment and the day before behavioral testing. DNA was isolated from a subset of subject stool samples for subsequent 16 S rRNA sequencing to quantify specific microbiota composition. A set of brain sections from all subjects in Experiment 1 were processed for immunohistochemical labeling of Iba-1, a microglia-specific marker, to quantify microgliosis in the PVN, NAcc, DG, and DR.

Results

Social isolation alters social affiliation and anxiety-like behaviors

In the SA test, no significant treatment effects were found in the total time spent in the conspecific cage (F1,40 = 0.45, p = 0.50; Table S1) nor in the total number of conspecific cage entries (F1,40 = 2.82, p = 0.10). There was a significant treatment effect for overall social behavior (i.e. grooming and huddling), where socially isolated voles showed a decreased duration of affiliative behavior toward the unfamiliar conspecific (F1,40 = 4.36, p < 0.05) compared to controls (Fig. 1A). However, socially isolated voles had significantly increased duration (F1,40 = 4.09, p < 0.05) and frequency (F1,40 = 4.36, p < 0.05) of sniffing the unfamiliar conspecific compared to cohoused controls (Fig. 1C–D). Sex differences were also found: male voles showed increased rearing (F1,40 = 6.79, p < 0.05) and effort biting (F1,40 = 6.28, p < 0.05) in comparison to females, whereas females showed increased duration of avoidant behaviors to the conspecific (F1,40 = 4.41, p < 0.05) and general sniffing behavior (F1,40 = 11.03, p < 0.01) compared to males (Table S1). No other main effects of sex nor sex-by-treatment interactions were found.

Fig. 1

Social affiliation and anxiety-like behaviors were assessed using the SA and EPM tests. Socially isolated prairie voles spent significantly less time affiliating with an unfamiliar, same-sex conspecific during the SA test (A). There were no differences in the frequency of social behavior (B). Social isolation increased the duration (C) and the frequency (D) of sniffing behavior toward the unfamiliar, same-sex conspecific. Male and female prairie voles that were isolated for 6 wks spent a significantly less percentage of time on the open arms of the EPM (E). There were no differences in overall locomotor activity on the EPM (F). As there were no sex by treatment interactions, data were pooled from males and females. Bars indicate mean ± SEM. * represents p < 0.05. ** represents p < 0.01.

In the EPM test, socially isolated male and female voles spent a significantly lower percentage of time in the open arms of the EPM (F1,40 = 9.16, p < 0.01) compared to cohoused controls (Fig. 1E), but no group differences were found in locomotor activity (Fig. 1F). No overall sex differences nor sex-by-treatment interactions were found in any measured behaviors in the EPM test.

Social isolation alters neuronal activation in a sex- and brain region-specific manner

Immunoreactive staining resulted in Egr-1 labeling in selected brain areas (Fig. 2, Fig. 3). A main treatment effect was found in the DR, where isolated voles had a lower number of Egr-1 cells compared to cohoused controls (F1,40 = 6.09, p < 0.05). A significant sex-by-treatment interaction was found in the NAcc, where isolated females had higher cell counts of Egr-1 compared to cohoused females as well as both isolated and cohoused males. No treatment effects were found in the DG (F1,39 = 0.05, p = 0.83), suggesting a brain region-specific effect of social isolation. We also found overall sex differences, where males had higher Egr-1 counts in the DG (F1,39 = 8.67, p < 0.05) and PVN (F1,39 = 4.94, p < 0.05), but lower Egr-1 counts in the DR (F1,40 = 4.40, p < 0.05) compared to females. No other sex-by-treatment interactions were found in any of the quantified brain regions.

Fig. 2

Social isolation alters Egr-1 staining in a sex- and brain region-specific manner. Isolated female prairie voles had significantly higher levels of Egr-1 in the NAcc. Social isolation was associated with lower Egr-1 levels in the DR. Females had higher Egr-1 levels in the DR, but lower Egr-1 levels in the DG and in the PVN compared to males. Egr-1, early growth response protein 1, NAcc, nucleus accumbens, DR, dorsal raphe, DG, dentate gyrus of the hippocampus. Bars indicate mean ± SEM. * represents a main effect of sex, p < 0.05. # represents a main effect of treatment, p < 0.05. Bars with different letters differ significantly from each other.

Fig. 3

Representative images of immunostaining of Egr-1 in the NAcc, DR, and DG. Treatment and sex are depicted above each image. NAcc, nucleus accumbens, AC, anterior commissure, DR, dorsal raphe, DG, dentate gyrus of the hippocampus. Scale bar = 100 μm.

Social isolation affects OT neuronal activation in the PVN of female, but not male voles

Given the previous literature showing social isolation-altered OT activity in the PVN (Grippo et al., 2007b) and the significant sex difference in Egr-1 expression in the PVN found in the present study, we examined how social isolation altered OT cells and their activation in the PVN (Fig. 4). A significant sex-by-treatment interaction was found for total OT cells (F1,39 = 5.69, p < 0.05), where cohoused females had lower numbers of OT cells compared to isolated females and both isolated and cohoused males (Fig. 4D). A significant sex-by-treatment interaction was also found in the percentage of Egr-1 cells co-labeled for OT (F1,39 = 8.33, p < 0.01), where isolated females had a higher percentage of Egr-1/OT double-labeled cells over total Egr-1 cells compared to all other groups (Fig. 4E). A similar trend was also found in Egr-1/OT over total OT cells in the PVN, but this difference did not reach statistical significance (Fig. 4F). Sex differences were also found in the number of AVP cells (Table S2) and CRF cells (Table S3) in the PVN, but not of 5-HT cells in the DR (Table S4).

Fig. 4

Representative images of immunostaining of Egr-1, OT, and double-labeled Egr-1/OT cells in the PVN (A-B). Social isolation alters activation of OT system in a sex-dependent manner. Females had significantly lower levels of Egr-1 in the PVN compared to males (C). Isolated (Iso) females had significantly higher numbers of OT cells compared to cohoused (CH) control females (D). Isolated females also had a significantly higher percentage of Egr-1 cells co-labeled for OT (E). There were no significant differences in the percentage of OT cells co-labeled for Egr-1 (F). OT, oxytocin, Egr-1, early growth response protein 1, PVN, paraventricular nucleus of the hypothalamus. Bars indicate mean ± SEM. Bars with different letters differ significantly from each other. * represents p < 0.05. Scale bar = 10 μm, 50 μm.

Social isolation alters neurochemical expression in a brain region-specific manner

To explore the neurochemical circuit involved in social isolation, we examined additional neurochemical marker expression by using Western Blotting in key brain regions known to be involved in social behaviors. Our data show that socially isolated voles had increased CRF in the NAcc (F1,38 = 4.68, p < 0.05) as well as increased BDNF (F1,31 = 4.06, p < 0.05) and OTR (F1,36 = 5.53, p < 0.05) in the AMY in comparison to cohoused controls (Fig. 5). An overall sex difference was also found, females (3.63 ± 0.74) had increased levels of CRFR2 in the PVN compared to males (1.80 ± 0.16) (F1,37 = 6.01, p < 0.05) (Table S5). No main effects of treatment, sex, or sex-by-treatment interactions were found in other selected neurochemical markers/brain regions (Table S5).

Fig. 5

Social isolation alters neurochemical markers in a brain region-specific manner. Isolated animals had increased levels of CRF in the NAcc (A) and BDNF and OTR in the AMY in comparison to cohoused control subjects (B). Neurochemical marker expression levels are plotted as over individual housing marker GAPDH to serve as a loading control. As no sex by treatment interactions were found, data were pooled from males and females. CRF, corticotropin releasing factor, NAcc, nucleus accumbens, BDNF, brain derived neurotropic factor, OTR, oxytocin receptor, AMY, amygdala, GAPDH, glyceraldehyde 3-phosphate dehydrogenase. Bars indicate mean ± SEM. * represents p < 0.05.

Social isolation alters microgliosis in a brain-region and sex-specific manner

Immunostaining for Iba-1 resulted in specific labeling in selected brain areas in prairie voles (Fig. 6A–H). A treatment effect was found in the DR, where socially isolated voles had significantly higher numbers of Iba-1 cells compared to cohoused controls (F1,40 = 4.38, p < 0.05). A sex difference was also found in the DR (F1,40 = 14.90, p < 0.01) and in the PVN (F1,42 = 69.51, p < 0.01), where males had higher Iba-1 cell counts compared to females. No sex-by-treatment interactions were found in these regions. Significant sex-by-treatment interactions were found in both the NAcc (F1,42 = 7.96, p < 0.01) and DG (F1,37 = 4.94, p < 0.05). Post-hoc analyses indicated that socially isolated females had increased Iba-1 cell counts in both the NAcc and DG compared to cohoused females; cohoused females also had significantly less Iba-1 cells compared to cohoused males in the NAcc (Fig. 6I).

Fig. 6

Representative images of immunostaining for Iba-1 in the NAcc, DR, PVN, and DG (A-H). Social isolation alters microgliosis in a sex- and brain region-dependent manner (I). Isolated females had significantly higher levels of microgliosis, as indicated by increased Iba-1 cell counts, in the NAcc and in the DG compared to cohoused control females. Isolated subjects had significantly higher levels of microgliosis in the DR. Males had significantly higher levels of microgliosis in the DR and in the PVN in comparison to females. Iba-1, ionized calcium binding adaptor molecule 1, NAcc, nucleus accumbens, DG, dentate gyrus, DR, dorsal raphe, PVN, paraventricular nucleus of the hypothalamus. Bars indicate mean ± SEM. Bars with different letters differ significantly from each other. ** represents a main effect of sex, p < 0.01. # represents a main effect of treatment, p < 0.05. Scale bar = 100 μm.

Social isolation alters CORT levels in a sex-specific manner

We processed plasma samples for CORT to examine how social isolation altered physiological stress responses. We found a significant sex-by-treatment interaction (F1,87 = 5.19, p < 0.05), where socially isolated males had higher CORT levels compared to male controls, but this effect was not seen in females (Fig. S1A). When looking at the data separated by experiments, we found the same sex-by-treatment effect in EPM animals (F1,42 = 5.27, p < 0.05), where isolated males had increased CORT compared to control males (Fig. S1B). Although not significant (F1,41 = 0.92 p = 0.34), the same pattern was found in SA subjects (Fig. S1C).

Social isolation alters gut microbial diversity in male and female prairie voles

Data from the 16s rDNA analysis (Fig. 7) indicated that the most abundant phylum was Bacteroidetes (56%), with the majority of the Bacteroidetes in the family Muribaculaceae/S24-7 (45%) and the family Prevotellaceae (8%). The Firmicutes phylum was next in abundance (33%), with most of this split among the Ruminococcaceae (15%), the Lachnospiraceae (12%), and the Erysipelotrichaceae (4%). The taxonomic composition of the M. ochrogaster microbiome observed in this study using the V3–V4 16 S rDNA amplicon marker is in overall agreement with the composition observed with the V4 region in a previous study (Curtis et al., 2018) that compared the taxonomic distribution measured with different 16 S amplicons.

Fig. 7

Mean relative abundance of bacterial taxa in stool from the different groups at the phylum (A) and family (B) level. Bacteroidetes was the most abundant phylum in all groups (56% overall, orange bars in (A)). Muribaculaceae (formerly S24-7) was the most abundant family in all groups (45% overall, orange bars in (B)). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

The only significant difference in alpha diversity between groups at the family level was between pre-vs post-cohousing samples in the InvSimpson's index (p = 0.031), which emphasizes the evenness of the distribution of the taxa over the richness of taxa. There were no significant differences in alpha diversity at the species level (see Table S6 for all alpha-diversity values).

Differences in beta-diversity for the DESeq2-normalized taxa counts from the samples were tested using the Bray-Curtis dissimilarity (Bray and Curtis, 1957), Unweighted Unifrac (Lozupone and Knight, 2005) and Weighted Unifrac (Lozupone et al., 2007) (see Table 1 for all beta-diversity values). PERMANOVA analysis with adonis2 (Oksanen et al., 2015) found significant differences in Bray-Curtis dissimilarity at the species level between isolated and cohoused post-treatment samples (p = 0.023), and for the pre-vs. post-treatment samples for cohoused animals (p = 0.033). There were no significant differences in beta diversity for pre- vs. post-treatment samples in the isolated animals. For Weighted Unifrac, social environment manipulation for post-treatment samples was significant at the species level (p = 0.034), but not the family level. In addition, pre- vs. post-treatment stool analyses for cohoused animals was significant at both the species (p = 0.041) and family levels (p = 0.031). Significance for Bray-Curtis and Weighted Unifrac, but not Unweighted Unifrac, suggests that the differences are in abundant taxa. An additional analysis of the samples in 3 groups (all pre-treatment vs post-isolation vs post-cohousing) was significant only for Weighted Unifrac at the species level (p = 0.024), suggesting that the differences are in abundant lineages that are phylogenetically diverged from one another. The contribution of variables to total variation between samples at the species level was also analyzed using CAP ordination (Fig. 8). For the post-treatment samples, social manipulation contributed significantly to variation determined by Weighted Unifrac (p = 0.022) (Fig. 8A). The analyses for isolated animals are shown in Fig. 8B, and analyses for the cohoused animals are shown in Fig. 8C. When the pre- and post-treatment samples were compared for each group separately, the model generated in the CAP analysis indicates that housing manipulation (pre-vs. post-treatment) did not contribute to the change for isolated animals (Fig. 8B) but did contribute to the change for cohoused animals by both Bray-Curtis (p = 0.001) and Weighted Unifrac (p = 0.008) (Fig. 8C). Thus, the CAP results and the PERMANOVA adonis2 results are in agreement that pair-housing contributes significantly to the beta diversity of post-treatment samples as well as to the post-treatment change in diversity for cohoused animals.

Table 1

PERMANOVA of the effects of variables on beta-diversity determined using adonis2.

Variables		Stool Samples		Bray-Curtis		Unweighted Unifrac			Weighted Unifrac
		R2			P	R2	P		R2		P
Species-level differences	Sex	All	0.011		0.386	0.014		0.171		0.015		0.157
	Social manipulation	all	0.015		0.131	0.01		0.541		0.011		0.454
	Pre- vs. Post-	all	0.012		0.339	0.008		0.758		0.013		0.262
	Social manipulation	pre-only	0.025		0.270	0.018		0.625		0.027		0.254
	Social manipulation	post-only	0.040		0.023*	0.015		0.820		0.049		0.034*
	Pre- vs. Post-	isolation	0.023		0.385	0.014		0.893		0.035		0.127
	Pre- vs. Post-	cohoused	0.035		0.033*	0.016		0.775		0.044		0.041*
	Pre-all vs Isolation-post vs Cohoused-post	all	0.031		0.064	0.017		0.826		0.039		0.024*
Family-level differences	Sex	all all allpre-only post-only isolation cohousedall	0.0140.010.0080.0310.0320.0330.0230.023		0.2570.4490.6200.2000.1900.1590.3330.330	0.0170.010.0080.0110.0060.0100.0060.019		0.1440.6520.6480.8340.9370.8210.9450.584		0.0180.0080.0110.0420.0450.0300.0670.035		0.1370.5430.4010.1210.1120.2380.031*0.128
	Social manipulation
	Pre- vs. Post-
	Social manipulation
	Social manipulation
	Pre- vs. Post-
	Pre- vs. Post-
	Pre-all vs Isolation-post vs
	Cohoused-post

* Significant at p < 0.05.

Fig. 8

Constrained Analysis of Principal coordinates (CAP) analysis. The ordination of the data from the beta-diversity metrics for (A) Post-housing-treatment samples, with cohoused animals shown in red and isolated animals shown in green; (B) Isolated animals only, with pre-housing-manipulation samples shown in green and post-housing-manipulation samples shown in red; (C) Cohoused animals only with the same data point color scheme as in (B). For all ordinations, the most parsimonious model formed is shown, with the p value of the model, and the proportion of inertia (variation) explained by constrained variables. CAP1 and CAP2 are the two greatest contributing components of the constrained variables. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Social isolation alters the abundance of specific microbial taxa in male and female prairie voles

The differences in abundance of taxa post-relative to pre-treatment were compared using DESeq2 (Table 2) where all species-level changes are shown that have a p < 0.1 adjusted for multiple comparisons. In isolated animals, an uncultured Anaeroplasma group decreased, whereas several other taxa, including Candidatus Saccharimonas uncultured rumen bacterium, Desulfovibrio uncultured bacterium, Lachnospiraceae UCG-001 metagenome species, Muribaculaceae and Gastranaerophilales gut metagenome species, increased post-isolation. Intriguingly, these same Candidatus Saccharimonas uncultured rumen bacterium and Lachnospiraceae UCG-001 metagenome species decreased from pre-to post-treatment in cohoused animals. Christensenellaceae uncultured species and Coprococcus 2 uncultured bacterium also showed decreases whereas a Clostridiales vadinBB60 group species showed an increase in cohoused animals. Various changes were also found in males and females under isolation and cohoused conditions (Table 2). All species-level whole-group changes with a log2fc |0.693| are shown in Tables S7–S12.

Table 2

Species that change by the whole-group method post-relative to pre-housing manipulation with p adj <0.1. Increases post-housing are shown in green, while decreases are shown in red. *top part of the cell is the log2 fold change post- relative to pre-housing manipulation. ^bottom part of the cell is the p-value adjusted for multiple comparisons.

In order to minimize the effect of individual differences prior to housing manipulation, differences in microbiota were also analyzed by determining the pre-to post-treatment change on a ‘vole-by-vole’ (VBV) pairwise basis (see methods). In isolated animals, a post-treatment reduction was found in an Anaeroplasma uncultured bacterium and a Candidatus Saccharimonas uncultured bacterium (Table 3). Isolation-induced reductions in these two along with additional taxa were found in isolated males and females, respectively. In cohoused animals, an increase in a Clostridiales vadinBB60 group uncultured bacterium was found in females and of an Alistipes uncultured bacterium was found in males (Table 3). All VBV fold change and log2 fold change values are shown for each taxon in Table S13.

Table 3

Significant differences at the species level between treatment groups on tested ‘vole-by-vole’. Increases post-housing are shown in green, while decreases are shown in red. All of the changes detected by the VBV method are in strains in the highest abundance quartile in pre-treatment animals. In contrast, the Anaeroplasma species was the only one detected by the whole-group method that is in the highest abundance quartile. *Top part of the cell is the number of animals out of total in which there is a pretreatment to posttreatment change in the same direction. ^Bottom part of the cell is the average per animal log2 (fold change) from pretreatment to posttreatment for all animals in the treatment group.

Microbial taxa correlate with neurochemical & microglial marker expression in the brain and behavioral phenotypes

Spearman Rho correlations were performed to examine significant associations across microbial abundance, neurochemical & microglial marker expression, and behavior (Table 4). As no significant correlations were found when all animals were grouped together, animals were then separated by treatment and sex. In isolated animals, significant correlations were found in both directions for several microbial taxa with OTR in the AMY or with Iba-1 in the DR. It is interesting to note that an increase in Prevotellaceae species negatively correlated with both OTR in the AMY and Iba-1 in the DR in isolated animals. In cohoused animals, significant correlations were found in both directions for several microbial taxa with CRF in the NAcc, OTR in the AMY, or with Iba-1 in the DG or PVN. There was virtually no overlap in significant correlations across the isolated and cohoused animals. Further, significant correlations in both directions between several microbial taxa and anxiety-like behavior or social affiliation were found only in the cohoused animals. Interestingly, an increase in an Allobaculum gut metagenome species was negatively correlated with open arm percentage, but positively correlated with social affiliation. Finally, several significant correlations were also found in males and females, respectively (Table 4).

Table 4

Microbial taxa correlate with neurochemical & microglial marker expression in the brain and behavioral phenotypes. Positive correlations are shown in green. Negative correlations are shown in red. *denotes rs values. ^ denotes p values. OA = open arm.

Discussion and conclusions

In the present study, chronic social isolation affected the gut-immune-brain axis and relevant behaviors in adult prairie voles. Given the highly social nature of prairie voles, the lack of social interactions is particularly detrimental; thus, social isolation is a well-established stress paradigm commonly used in prairie vole studies. Chronic social isolation has been shown to enhance HPA activity (as indicated by elevated levels of CORT, ACTH, and adrenal weight) in prairie voles, which provides substantial validation for using the chronic isolation paradigm to induce stress (Grippo et al., 2007a, 2007b; Jarcho et al., 2019; Ruscio et al., 2007; Watanasriyakul et al., 2019). Our data indicate that social isolation stress reliably increased anxiety-like behaviors in both males and females, which is consistent with previous prairie vole studies (Grippo et al, 2008, 2014; Lieberwirth et al., 2012) as well as in other rodent species (Harvey et al., 2019; Ieraci et al., 2016) and in humans (Domènech-Abella et al., 2019). Previous literature also shows that social isolation can either increase (Latane et al., 1972; Lieberwirth et al., 2012) or decrease (Okada et al., 2015; Shoji and Mizoguchi, 2011) social behaviors in various rodent models depending on species specificity and/or variations in experimental paradigms. Our data indicate that chronic isolation for six weeks led to a decrease in social affiliative behaviors, including social grooming and huddling, in adult prairie voles. As social stressors can increase both anxiety and social avoidance (Hammels et al., 2015; Iñiguez et al., 2014), decreased affiliation in our study may be induced either directly from social isolation, indirectly from increased anxiety-like and/or social avoidance behaviors, or from a combination of both. Interestingly, isolated prairie voles increased social sniffing toward unfamiliar conspecifics during the SA test. Social isolation can impair social recognition memory in other rodent species (Leser and Wagner, 2015; Liu et al., 2018; Shahar-Gold et al., 2013), increase social sniffing behavior (Shimozuru et al., 2008; Shoji and Mizoguchi, 2011), and alter olfaction independent of social memory (Gusmão et al., 2012). It would be interesting to further investigate our sniffing data by assessing whether social isolation impairs social recognition memory, social novelty/investigation, or both. Together, our data indicate that chronic isolation had detrimental effects on behaviors in male and female prairie voles.

Our data also indicate that chronic isolation altered neurochemical expression and neuronal activation in multiple brain areas. Social isolation increased CRF expression in the NAcc as well as BDNF and OTR expression in the AMY. These data are generally consistent with previous studies showing that social isolation alters CRF and OT systems in prairie voles (Pan et al., 2009; Pournajafi-Nazarloo et al., 2011; Ruscio et al., 2007), mice (Senst et al., 2016), and rats (Harvey et al., 2019; Oliveira et al., 2019) and can also affect the BDNF system in mice (O'Keefe et al., 2014). As CRF (Cipriano et al., 2016; Crumeyrolle-Arias et al., 2014; Kasahara et al., 2011; Pomrenze et al., 2019; Takahashi, 2001), OT (Cavanaugh et al., 2018; Gottschalk and Domschke, 2018; Li et al., 2016; Popik and Van Ree, 1999; Raam et al., 2017; Smith et al., 2016), and BDNF (Bahi, 2017; Berry et al., 2015; Miao et al., 2018; Scattoni et al., 2013; Zheng et al., 2016) systems have all been implicated in anxiety-like behavior, social affiliation, and social sniffing/recognition, alterations in those neurochemical systems may contribute to behavioral changes seen from social isolation in prairie voles. The DR has also been shown to be impacted by social isolation in mice (Lukkes et al., 2009; Sargin et al., 2016). Our data illustrates decreased Egr-1 labeling in the DR in socially isolated voles. As no effects were found in Egr-1/5-HT double labeling in the DR, this decreased neuronal activation might be occurring in non-serotonergic neurons. Interestingly, social isolation can alter dopaminergic neurons in the DR, and these changes mediate the alterations in social behaviors from isolation (Matthews et al., 2016). Collectively, our data support the notion that a neurocircuit of multiple brain areas and neurochemical systems is involved in mediating impacts of social isolation on complex behaviors (Cacioppo et al., 2015; Mumtaz et al., 2018).

Our sex-specific effect of social isolation on OT in the PVN should be noted. Socially isolated female voles had enhanced OT labeling and increased Egr-1/OT double labeling in the PVN compared to cohoused females, but this effect was not found in males. These data support the previous finding showing PVN OT alterations in isolated female, but not male voles (Grippo et al., 2007b) as well as female-specific alterations in PVN signaling in mice (Senst et al., 2016). Our data is also supported by previous findings in female prairie voles that increased anxiety-like behaviors are associated with increased OT expression/activity in the PVN (Smith et al., 2016; Smith and Wang, 2014). This sex- and brain region-specific effect of social isolation is further illustrated by increased NAcc Egr-1 expression in female, but not male, voles. It should be noted that the vole's NAcc receives OT projections from the PVN (Ross et al., 2009), is enriched with OTR (Insel and Shapiro, 1992), and serves as a critical node for social behavior (Tabbaa et al., 2017b).

The microgliosis data in our study are exciting – social isolation increased Iba-1 expression in the DR of both male and female voles as well as in the NAcc and DG in female voles, indicating a brain region- and sex-specific effect. These data offer a unique and novel look at how central glial cells may be shaping social circuits in the brain and resulting behaviors. Interestingly, the impacts of social isolation on Iba-1 and Egr-1 were in the same direction in the NAcc, but in the opposite direction in the DR. As microglia have bidirectional communication with neurons and can alter neuronal signaling (Eyo and Wu, 2013), our data suggest that glial cells may play differential roles in a brain region-specific manner. The causal role of microglia inducing behavioral and neuronal alterations from social isolation has been previously demonstrated in rats (Wang et al., 2017). Microglia also play a critical role in dopamine receptor elimination in the NAcc to shape sex-specific social behaviors in rats (Kopec et al., 2018).

When looking at our microbiome data, the most striking change in a specific taxon after social manipulation is the post-isolation reduction of an Anaeroplasma species detected by both the whole-group and the VBV methods. Though the log2fc is not very large by either method, it is a change detected by both methods in one of the most abundant taxa in the vole microbiome. Anaeroplasma are anaerobic, gut-dwelling mycoplasmas (Joblin and Naylor, 2002) that have recently been shown to correlate with high levels of IgA+ plasma cells in the intestinal mucosa of mice (Beller et al., 2020). Transfer of Anaeroplasma-enriched filtrate from fecal material of high IgA BALB/c mice to antibiotic-treated, microbiota-ablated mice was sufficient to induce an increase in fecal IgA levels, indicating a beneficial role for Anaeroplasma in mucosal immunity. Stress-induced social avoidance is correlated with changes in gut microbiome composition and immune function in mice as well (Szyszkowicz et al., 2017). The reduction of microbiota associated with proper immune function along with the brain region-specific alterations in microgliosis in isolated prairie voles emphasize the potential role of the gut-immune-brain axis in mediating isolation-induced outcomes.

By the whole-group method, a Ruminococcaceae UCG-014 species decreased in isolated females whereas a Butyrivibrio species decreased in isolated males, indicating a potential sex-specific role. Ruminococcaceae UCG-014 has been reported to be at a lower abundance in humans with higher anxiety (Chen et al., 2019). Strains of Butyrivibrio are producers of the short-chain fatty acid butyrate, which can suppress inflammatory responses (Forbes et al., 2016; Singh et al., 2014). Butyrivibrio is also associated with higher quality of life scores in a large study of microbiota and depression in humans (Valles-Colomer et al., 2019). Further, it is intriguing to note that isolation treatment also led to increases in some of microbial taxa such as Lachnospiraceae UCG-001 and Desulfovibrio. Species of these genera have been found enriched in neonatal rats after maternal separation (Rincel et al., 2019), suggesting a possible role for these species in negative consequences that arise from reduced social engagement.

Global differences in microbiome diversity after housing treatment were apparent in cohoused, but not isolated, animals. Our data show that 1) physiological changes due to isolation did not significantly alter beta-diversity and 2) physiological effects of cohousing and/or sociality in a novel environment led to significant beta-diversity differences between the baseline and final stool collection (pre- vs. post-cohousing samples). The significant difference in beta-diversity in the cohoused animals before and after the cohousing period could be due to physiological effects of cohousing over time, the transmission of microbes from cagemates while in a novel environment, or a combination of both. It should be noted that there were no novel social interactions in the control group, as the cohoused subjects were already housed with their cagemates prior to the start of the experiment. Therefore, enriched microbiome diversity in cohoused, but not isolated, animals may indicate that sociality in a novel environment may play a role in modulating microbiome diversity. Perhaps the presence of a familiar conspecific and resulting social interactions in novel environments lead to fluctuations in the gut microbiome that may be adaptive, whereas the lack thereof does not – it would be interesting to further explore these ideas in subsequent studies.

We performed correlations to further explore relationships between microbial alterations with the brain and behaviors. Several significant correlations are worth noting. For example, an increase in an Oscillibacter 1–3 species was positively correlated with OTR in the AMY in isolated animals. The abundance of this microbial species increased in males but decreased in females after social isolation, suggesting a complex interaction between the isolation-induced stress, sex and the abundance of this species. The roles Oscillibacter plays in gut health are not yet fully understood, although this genus has been observed to increase in abundance in high-fat diet-induced obesity in mice and to positively correlate with gut barrier disfunction (Lam et al., 2012). Also, the same Gastranaerophilales bacterium that decreased in abundance after social isolation in males also had a strong negative correlation with OTR in the AMY in females. Since OTR expression in the AMY was increased in isolated animals, this result, combined with the differential abundance results, suggest that this Gastranaerophilales species might be associated with lower levels of stress. Interestingly, this group is also more abundant in mice with reduced inflammation phenotypes after chronic unpredictable stress relative to controls, suggesting an association of Gastranaerophilales with resilience in response to stress (Inserra et al., 2019). A negative correlation between an increase in Prevotellaceae UCG-001 bacterium P3 with Iba-1 in the DG and the PVN in cohoused animals is also worth mentioning. This was somewhat surprising as Prevotellaceae UCG-001 is one of the genera with higher abundance in isolated rats (Dunphy-Doherty et al., 2018), but microgliosis was not measured and this effect might also be species-specific. Another study showed that prebiotic treatment can result in Prevotellaceae UCG-001 enrichment, leading to increased short-chain fatty acids (SCFAs) (Song et al., 2019). Given the ability of SCFAs to interact with the immune system, a further exploration of the link between Prevotellaceae UCG-001 and central immune modulations is worth noting. We also found a positive correlation in cohoused animals between social behavior and an increase in Ruminococcaceae UCG-010 uncultured bacterium, suggesting a potential microbial target for driving prosocial behavior. Interestingly, we found that an increase in an Allobaculum gut metagenome species was negatively correlated with open arm time, but positively correlated with social affiliation. Further, there were many other microbial species found to correlate with neurochemical expression, microgliosis marker, and behaviors in a marker-, brain region-, sex-, and behavior-specific manner with limited overlap. These data may suggest a role of dynamic interactions and relative abundances of the microbial community in modulating complex bodily functions and resulting behaviors (Lin and Zhang, 2017; Mayer et al., 2015; Rooks and Garrett, 2016).

The inclusion of both male and female subjects in the current study revealed remarkably interesting sex differences and sex-specific brain and gut microbial alterations in response to social isolation. For example, in females, social isolation altered Egr-1 and Iba-1 labeling in the NAcc – a brain region in which microglia have been shown to shape neurochemical circuits underlying sexually dimorphic social behaviors (Kopec et al., 2018). Increased Egr-1/Iba-1 labeling were also found in the PVN in females, although their causal relationship still needs to be tested. As mentioned, the role of OT on anxiety-like and social behaviors have been well documented; our data illustrate both female-specific PVN OT activation in response to isolation as well as a positive correlation between an increase in a Bacteroidia species and OTR expression in the AMY, indicating a potential role of gut microbial changes in differential oxytocinergic system activation in females. Although a majority of the sexually dimorphic isolation effects found in the brain were in females, we also found that social isolation significantly elevated CORT in male, but not female voles. Interestingly, previous studies have shown sexually dimorphic effects of CORT on social behavior in prairie voles (Devries et al., 1996). Further, our data also show male-specific changes in microbiome alterations. For example, certain microbial taxa, such as Oscillibacter sp. 1–3 species, increased in isolated males, but decreased in isolated females. Previous research have consistently shown bidirectional crosstalk between gut microbiota and the HPA axis (Farzi et al., 2018). In one study, male rats that underwent prenatal stress had increased Oscillibacter, and changes in the microbiome were correlated to changes in stress responsivity (Golubeva et al., 2015). In another study, probiotic administration was sufficient in lowering stress-induced CORT in male mice, which supports the notion that microbial changes can alter HPA activity from stress (Bravo et al., 2011). Further, housing density stress increased CORT in males of another vole species, and CORT levels were correlated with gut microbial alterations as well (Liu et al., 2020). Overactivation of the HPA axis has been linked to increased anxiety-like behavior, altered neuroimmune responses, and altered microbiota composition, indicating that these systems commonly interact to shape behavior (Amini-Khoei et al., 2019). These drastic differences in microbial, neurochemical, and physiological data across sexes in our study are intriguing and further emphasize the need to include both male and female subjects in subsequent causal explorations of these underlying mechanisms of social isolation.

Finally, it is worth noting that we did not see sex-specific alterations in the tested behaviors, which is consistent with previous prairie vole studies (Grippo et al., 2007b; McNeal et al., 2014). It is well recognized that sexual dimorphisms in physiology and neural substrates may underlie sex differences in behaviors, including stress responses (Heck and Handa, 2019; Nelson and Lenz, 2017; VanRyzin et al., 2018). It is also recognized that sexually dimorphic neurochemical systems may allow males and females to have compensatory mechanisms that work in concert with their physiology to produce similar behavioral outcomes (De Vries and Forger, 2015; De Vries and Villalba, 1997; Grabowska, 2017; Ross et al., 2009). Interestingly, our data show that social isolation altered Egr-1, OT, and Iba-1 in selected brain areas in females, but altered circulating CORT in males. Therefore, social isolation stimuli may affect anxiety-like and social behaviors via brain and endocrine systems differentially between males and females, and specific microbial taxa may be involved in modulating these processes, as gut microbiota have been shown to be involved in shaping sexual dimorphisms in physiological systems (Jaggar et al., 2020; Jašarević et al., 2016; Thion et al., 2018). Taken together, our data suggest that social isolation can alter physiological and neurochemical systems in a sex-dependent manner. As isolation altered central measures in females, yet physiological measures in males, perhaps future studies should differentially assess pharmacological targets in male and female subjects to better understand the functional implications of these changes. This interesting mismatch between behaviors, neural substrates, microgliosis, and microbial alterations in the present study indicate both an opportunity and a necessity to further explore the role of differential activation of the gut-immune-brain axis and neuronal circuits on the various negative outcomes seen from social isolation in male and female prairie voles. Given the variety of significant changes found in the current study, there should be future investigations which assess both peripheral and central mechanisms and how they are causally related to one another. Future studies may utilize pharmacological techniques, such as the use of minocycline to block region-specific activation of microglia and/or neurochemical systems, or targeted probiotic intervention to directly boost reduced microbial populations, such as Anaeroplasma, to better understand the functional implications of the changes seen in the current study.

CRediT authorship contribution statement

Meghan Donovan: Conceptualization, Investigation, Formal analysis, Visualization, Writing - original draft, Writing - review & editing, formation, review & editing. Calvin S. Mackey: Investigation, Formal analysis, Visualization, Writing - review & editing. Grayson N. Platt: Investigation. Jacob Rounds: Investigation. Amber N. Brown: Investigation. Darryl J. Trickey: Investigation. Yan Liu: Investigation, Formal analysis. Kathryn M. Jones: Visualization, Data curation, Supervision, Validation, Writing - original draft, Writing - review & editing, formation, review & editing, Funding acquisition, Resources. Zuoxin Wang: Conceptualization, Visualization, Supervision, Methodology, Data curation, Validation, Writing - original draft, Writing - review & editing, formation, review & editing, Funding acquisition, Project administration, Resources.

Declaration of competing interest

All authors have no declared conflict of interest and have nothing to disclose.