Behavior
Detailed statistical analysis is presented in Table S1.
Nesting All mice exhibited similar nesting behavior, with relatively high nesting scores regardless of Genotype, Group or Sex. No main effects or interactions were found (Fig S1).
Open fieldMales and females differed in total locomotor activity (Sex main effect; males>females) and in the time spent in center over time (Time*Sex interaction). Therefore, male and female data were analyzed separately. In male mice, all groups exhibited less locomotor activity over time, and C-Glud1+/- mice travelled less distance compared to C-Cre+ controls, regardless of stress exposure, during the first 40 minutes (Fig 1aLeft) and over the entire test period (1aRight). No Group effect or Genotype*Group interaction were found. In female mice, all groups exhibited less locomotor activity over time (1bLeft). We found no effect for Genotype, Group or a Genotype*Group interaction (1bRight).
Examination of time spent in the center of the open field in males revealed a main effect of Time and a Time*Genotype interaction. We found a significant change over time for C-Glud1+/- mice, which exhibited linear and polynomial trends, indicative of ‘erratic’ behavior (1cLeft). There was no change over time in C-Cre+ mice. During the first 5 min of the test (bin 1, 1cRight), C-Glud1+/- mice spent less time in the center than their C-Cre+ controls, and stress-exposed mice spent less time in the center than their stress-naïve controls. No Genotype*Group interaction was found. In females, we found a main effect for Time, but no Time*Genotype interaction (1dLeft). Similarly to male mice, C-Glud1+/- females spent less time in the center compared to their C-Cre+ littermates, and stress exposure similarly led to a decrease in time in the center (1dRight). No Genotype*Group interaction was found.
Social behaviors
We found no Genotype, Group or Sex effects, nor interactions between them, in Social Preference or Social Recognition (Fig S1).
Water T-maze
Male and female data were combined since no Sex effects or interactions were found.
Acquisition phase –Genotype*Group interactions for both Trials to Criterion and Reaction Time variability (RTv) were found. Post-hoc tests for both measures indicate that C-Glud1+/-/ Stress mice required more trials to criterion (1eLeft) and displayed higher RTv (1eright; Fig S2) than mice in all other groups. Data from 1 C-Glud1+/-/Stress mouse was excluded (>2SD above mean) in the RTv analysis. RT data followed the same pattern of results: RTs were highest in the C-Glud1+/-/ Stress group (not shown, see Table S1). Four mice from the C-Glud1+/-/Stress group [21.05%; significantly more than in other groups (0%)] failed to achieve criterion in the acquisition phase and received a maximal score of 60. These mice did not continue to the probe or reversal phase. Reversal phase – As in the acquisition phase, there was a Genotype*Group interaction for both Trials to Criterion and RTv. Post-hoc tests showed that C-Glud1+/-/ Stress mice required more trials to criterion (1fLeft) and displayed higher RTv (1fRight) than mice in all other groups. Higher RTs were observed in this group as well.
RNA-seq
Differentially expressed genes (DEGs) – as can be seen in Table 1, RNA-seq analysis showed that 5,641 genes (36% of the protein-coding genes) in the C-Glud1+/-/ Stress group, 88 in the C-Glud1+/-/ Control group, and 4 in the C-Cre+/ Stress group were differentially expressed compared to C-Cre+/Control (FDR-corrected q value <.10).
Table 1: Summary of protein-coding genes detected, tested and differentially expressed in the 3 experimental conditions. Tested = genes analyzed for differential expression (i.e., had normalized counts above 10 for at least one sample in each group).
Comparison (vs C-Cre+/Control)
|
Total Number of Detected Genes (protein coding only)
|
Tested
|
Significantly Upregulated
|
Significantly Downregulated
|
C-Cre+/ Stress
|
20,253
|
15,708
|
0
|
4
|
C-Glud1+/-/ Control
|
20,253
|
15,705
|
13
|
75
|
C-Glud1+/-/ Stress
|
20,253
|
15,640
|
2,837
|
2,804
|
As can be seen in Fig 2, in the C-Cre+/ Stress group all 4 DEGs were downregulated (2a), in the C-Glud1+/-/ Stress group the majority of the 88 DEGs were downregulated (2b), and the Glud1+/-/Stress group exhibited the highest number of DEGs, with an approximately equal number of upregulated and downregulated genes (2c). Examining the overlap between DEGs in our groups (2d), we found that all DEGs in the C-Cre+/ Stress group and most DEGs in the Glud1+/-/ Control groups were also downregulated in the Glud1+/-/ Stress group. There was only a small overlap between the C-Glud1+/-/ Stress group and other groups: out of 5641 DEGs, 5568 (98.7%) were unique to the C-Glud1+/-/ Stress group.
To detect enriched biological categories, we performed an over-representation analysis (ORA) on the 5641 DEGs in the C-Glud1+/-/Stress group using Enrichr46,47. As can be seen in Fig 3a, a BioPlanet 2019 ORA revealed significant changes in glutamate- and GABA-related pathways. These pathways were also prominent when other gene databases were used (e.g., GO Biological Process, Elsevier Pathway Collection; Fig S3a-b). The Elsevier Pathway Collection also pointed to the CRH secretion regulation pathway.
We then proceeded to examine expression patterns of genes from over-represented categories in all experimental groups. C-Glud1+/-/Stress mice consistently showed expression changes in glutamate- (3b) and GABA-related genes (3c). We further constructed a list of functionally significant genes that contribute to homeostasis at the glutamate tripartite synapse, the E/I balance and cognitive function8,48, using the GO terms "Glutamatergic synapse" (GO:0098978) and "GABA-ergic synapse" (GO:0098982). Interestingly, Glud1 expression was not significantly altered, but many of the genes that contribute to glutamate and GABA homeostasis were differentially expressed in the C-Glud1+/-/Stress group, but not in C-Cre+/Stress or C-Glud1+/-/ Control mice (S4a-b). Three genes associated with HPA axis function were also differentially expressed exclusively in this group (S4c).
The STRING database49 was used to examine protein-protein interactions between glutamatergic, GABAergic and stress-related DEGs (Fig S4). This analysis revealed a functional network with 7 clusters: Glutamatergic signaling (Red; 23 genes), G-protein transmembrane signaling (Brown; 11 genes), GABAergic (Dark green; 9 genes), neurotransmitter synthesis and packaging (light Green; 8 genes), CREB1 phosphorylation through NMDA receptor-mediated activation of RAS signaling (Teal; 4 genes), Potassium inwardly-rectifying channels (Blue; 3 genes) and RAS-ERK signalling pathway (Purple; 2 genes). As can be seen in Fig 3d, proteins encoded by DEGs interact within and between clusters.
We also examined whether genes known to interact with Glud1 (STRING) were differentially expressed in C-Glud1 +/- stressed mice. Out of 63 genes known to interact with Glud1, 25 were differentially expressed in C-Glud1+/-/Stress mice (16 upregulated).Protein-protein interactions between these genes are shown in Fig S5.
To assess the relationship between mPFC gene expression and behavior, we calculated Spearman correlations between DEG transcript levels and behavioral measures in the open field and T-maze tasks, in all mice. Thirty eight glutamatergic, GABAergic and stress-related genes differentially expressed exclusively in the C-Glud1+/-/Stress group were selected (3b-d, S4). As can be seen in Table S3, transcript levels correlated with several behavioral measures, but mostly with individual RTv in the Water T-Maze Reversal task. RTv values predicted transcript levels of 33 out of the 38 genes. Fig 4a depicts select correlations between glutamatergic, GABAergic and stress-related DEGs and individual RTv (higher RTv reflects compromised cognitive performance).
K-means cluster analysis on these genes (jamovi (v2.2.5) 'snowCluster' module (v6.7.0)), excluding 3 genes due to low factor loading (Standardized;<|.40|) with either cluster), yielded two clusters (4b). The final model (Table S4) accounted for 71.3% of the variance, with the first PCA accounting for 59.1% of the variance. The 2 clusters separated the samples into two non-overlapping groups (4c), which differ in expression levels of clustered genes (samples in the ‘blue’ cluster have high expression levels of genes marked in blue and low expression levels genes marked in orange, and vice versa).Notably,theseclusters differ in Reversal RTv scores (Mann–Whitney test; Z=-3.36, p<.001; 4d, Table S2). High- and low- performance mice (HP and LP, respectively) displayed differential expression of 34 of the 35 genes in the analysis; the top 12 differentiating genes are presented in 4e (Mann–Whitney tests; all Zs>3.31, all p’s<.001). Notably, the majority of LP/ Cluster 2 mice were C-Glud1+/-/Stress mice (orange; 60%; 33.33% of total sample), and the most differentiating gene was the mitochondrial glutamate carrier 1 gene, Slc25a22.
miRNA expression
We first conducted an in-silico search (using the DIANA miRNA database) for targets of our selected miRNA molecules (shown in the literature to be affected by stress and to be associated with cognitive function) which overlap with DEGs that emerged from the RNA-seq data. We found an approximately 40% average overlap between DEGs and targets of each miRNA. As can be seen in Fig 5a, several ofmiRNA targets that were also differentially expressed in C-Glud1+/-/Stress mice were shared by the selected miRNAs. The miRNAs examined were found to target genes related to glutamate, GABA and postsynaptic components in pathway analysis conducted with GO biological and cellular databases; several of the miRNAs were found to affect the same pathways (5b).
Next, we asked whether the expression of these miRNAs is altered in the mPFC of C-Glud1+/-/Stress mice. Detailed statistical findings are provided in Table S2. We examined 2 time points in separate cohorts: immediately following stress exposure (T1) and approximately 1 month later (T2; the timepoint at which the RNA-seq experiment was performed). Notably, RT-PCR analysis revealed that several of the DEGs that emerged at T2 were unaltered in the mPFC at T1 (Fig S6). As can be seen in Fig 5c, the expression of MiR203-5p increased by approximately 13-fold in mPFC at T1 in C-Glud1+/-/ Stress mice and remained mildly increased at T2 relative to both C-Cre+/Stress and C-Cre+/Control (dotted line) groups. The expression of MiR493-5p differed significantly between groups at T1, with higher expression in C-Cre+/ Stress mice and marginally higher expression in the C-Glud1+/-/ Stress group; these differences did not persist at T2.
Differences in the expression of MiR137-5p, MiR203-3p and MiR34a-5p emerged at T2. MiR137-5p expression was upregulated in the C-Glud1+/-/ Stress group compared to C-Cre+/ Control mice, and marginally compared to the C-Cre+/ Stress mice. MiR203-3p and MiR34a-5p were downregulated in C-Glud1+/-/ Stress mice compared to C-Cre+/ Control and C-Cre+/ Stress mice. For all 3 miRNA molecules, no differences were found at T1. No differences in MiR137-3p or MiR493-3p were detected at either time point. The mean miRNA expression score (mean ∆∆Ct across miRNAs for each mouse) at T1 was higher in C-Glud1+/-/Stress compared to the C-Cre+/ Stress (fold change: 2.60±0.24) and C-Cre+/ Stress (fold change: 1.21±0.49) mice. No such difference was found at T2.
Analysis of individual miRNA molecules and the mean miRNA expression score at T2 in the dCA1 and vCA1 hippocampal subregions revealed no changes in the C-Glud1+/-/ Stress group. In dCA1, the mean miRNA expression score and the expression of MiR137-5p were upregulated in the C-Cre+/Stress group (Fig S7).
Since mPFC MiR203-5p expression increased dramatically at T1, we asked which pathways are associated with targeted transcripts of this miRNA molecule. Pathway analysis revealed targets associated with cortisol synthesis, steroid metabolism, glutamate homeostasis and synaptic function (5d).
Finally, we examined the correlations between miRNA expression and behavioral measures in the water T-maze (5e; see additional behavioral measures in Table S3). MiR203-5p expression was correlated with the number of trials to criterion at both phases, but not with variability measures. In contrast,MiR34a-5p was correlated with variability measures in both phases, and with trials to criterion in reversal. MiR493-5p was not correlated with any of the measures. Overall, T-maze reversal performance was most highly correlated with GABA, glutamate and stress-related gene expression, as well as with miRNA expression levels, compared to other behavioral measures (Table S3).