Twenty-eight male mice (9 weeks old) of 28 different strains were purchased from Jackson Laboratory (Bar Harbor, ME, USA). Animals included in this study represent only a single cohort of mice that were randomly selected from a larger behavioral genetics study (Carhuatanta et al., 2014; Shea et al., 2015) using BXD recombinant inbred strains (For all strain information, please see Figure 2A), derived by crossing C57BL/6K (B6) and DBA/2J (D2), that have been significantly used to study gene-phenotype interactions. All animals included in this cohort arrived to our facility at the same time and completed the experimental procedures, including stressors, behavioral testing and sample collections, on the same day. All mice were singly housed with ad libitum access to food and water and allowed to acclimate to the Wright-Paterson Air Force Base (WPAFB) animal facility for 10 days. Ambient housing conditions were controlled for temperature (18–24°C), humidity (30–70%), and dark/light cycle (12:12; 06:00 on). All experiments were performed during the light phase. After completing the behavioral tests, mice were euthanized via rapid decapitation and whole brains were quickly frozen and stored at −80°C until further processing. All rodent handling and procedures were performed in accordance with the National Institute of Health standards and the Guide for the Care and Use of Laboratory Animals (National Research Council), and the protocol was approved by the WPAFB Institutional Animal Care and Use Committee. An overall research outline is shown in Figure 1.

Only stress groups were selected for this study because the purpose of this study was to investigate how cognitive performance under chronically unpredictable stress environments were various and how this variety was related to transcriptomic modification in the hippocampus. The chronic variable stress (CVS) paradigm was completed over 4 consecutive weeks as described previously (Schmidt et al., 2008; Carhuatanta et al., 2014). Briefly, mice were randomly exposed to twice daily stressors (morning and afternoon) for 5 days/week. Stressors included cold exposure (15 min at 4°C), hypoxia (30 min at 8–12% O₂), constant motion (60 min at 100 rpm), novel overnight housing, restraint (60 min) and open field (30 min in a 19″ × 10.5″ × 8″ cage) exposure. To limit predictability and habituation, stressors were presented in a randomized fashion to the mice.

All mice underwent Morris water maze (MWM) testing during the last week as described previously (Morris, 1984; Vorhees and Williams, 2006) to examine spatial memory acquisition. Briefly, a round basin (90 cm diameter) was used for MWM task and filled with water (42 cm of depth), mixed with nontoxic white paint until a clear platform became invisible (the range of water temperature: 19–24°C). The platform was located in the southwest quadrant submerged approximately 1 cm below the water surface and all visual distal clues were provided on the walls. Rodents were given 5 training days of 4 trials per a day using a random starting location. If rodents did not reach the platform within 60 s they were gently guided to the platform. The EthoVision video tracking software (version XT7.0.418, Noldus Information Technology, Leesburg, VT) was used to record mouse swim path, position, speed, and latency to platform. For this study, the mean of MWM performance (latency time to locate platform) during the 5 days was used to group the mice (Figure 2A) into high (better cognitive performers, shorter latency to platform time: < 20 s; n = 9) and low (worse cognitive performers, longer latency to platform time: >25 s; n = 9) MWM performance groups. Probe trial, for which the platform was removed, was also given to the animals on the 6th MWM testing day. Animals that displayed average MWM performance (latency to platform time: 20–25 s; n = 10) were excluded from the comparative analysis between high and low performing groups. The total distance moved during the MWM testing days and the probe trial was also recorded.

As soon as behavioral tests were completed (on the last day of 4 consecutive weeks for CVS), animals were immediately euthanized in about 10 min on the same day and, then, tissues were dissected and collected on ice and frozen quickly. Total RNA from hippocampus of the rodents were utilized for building library constructions and, then, running the Ilumina whole transcriptome RNA sequencer. The frozen brain was placed in the RNAlater®-ICE Frozen Tissue Transition Solution (Thermo Fisher Scientific, Waltham, MA) at −20°C for overnight and dissected hippocampus, and total RNA was extracted from the whole hippocampus using the RNeasy Mini Kit (Qiagen, Valencia, CA) according to instructions provided in the kit. RNA samples were shipped overnight to the UCLA Biosequencing core where the Illumina RNA sequencer was used for whole transcriptome RNA-sequencing. Samples were multiplexed, tagged with standard Illumina tags and sequenced using Illumina next-gen RNA-Sequencing by trained technicians.

All multiplexed data from the Illumina RNA Sequencer were downloaded and fastQ files of each sample were demultiplexed and merged into one fastQ file. All sequencing adapters were removed by running the FASTX-Toolkit, FASTQ clipper, and all reads were, then, trimmed by using another FASTX-Toolkit FASTQ trimmer (http://hannonlab.cshl.edu/fastx_toolkit/) if needed. The FastQC application (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to check quality control on FastQ files. All reads were aligned to UCSC genome browser (Mus musculus NCBIM37) by using TopHat (version 2.0.6) (Trapnell et al., 2012). The aligned reads (mean: 27.5 M reads) kept by TopHat running were quantified by HTSeq (version 0.5.3p9) (Anders et al., 2015). The DESeq package (Anders and Huber, 2010) in R was used to normalize read counts. Genes that had no reads across all samples were discarded by the DESeq analysis. The PANTHER (Protein ANalysis THrough Evolutionary Relationships) Classification System was used to classify genes in addition to use of the weighted gene coexpression network analyses (WGCNA) as described below.

PANTHER Classification System (http://PANTHERdb.org) (Thomas et al., 2003) was used to classify genes that were significantly expressed detected by the DESeq analysis and correlation analysis. Lists of Ensembl gene IDs were entered into the PANTHER database to classify genes into functional groups, providing information about their biological function and pathways. The PANTHER statistical over- and under-representation tests were carried out with log₂ fold change (log₂FC) values and available Entrez Gene IDs to investigate whether identified biological processes were statistically over- or under-expressed.

Normalized gene expression values from DESeq data were used to construct signed coexpression networks using the WGCNA package in R (Langfelder and Horvath, 2008). After the data cleaning step, the step-by-step network construction and module detection was carried out. The adjacency was calculated with the softPower of 10 (the scale-free topology fit index SFT R² = 0.815). The corresponding dissimilarity was calculated from topological overlap matrix (TOM) that was transformed from the adjacency, and then hierarchical clustering and dynamic tress cut methods were carried out to identify the modules with the recommended minimum module size of 30. The Dynamic Tree Cut method identified similar modules, their eigengenes were calculated, clustered based on their correlation and, then, merged modules if their correlation was > 0.75 (from 12 dynamic modules to 10 merged modules). An association of genes in each module with the traits (individual rodents' means of MWM performance during 5 days) was quantified by defining Gene Significance (GS), and a quantitative measure of module membership that is a correlation of the module eigengene and the gene expression profile was calculated. The WGCNA method identified modules with a statistically significant correlation between means of MWM performance and genes. The eigengene dendrogram function of the WGCNA was used to identify groups of correlated eigengenes, discovering mutual correlations between modules and the trait. Significantly correlated modules were used to perform GO enrichment analysis and identified the 20 best terms for each module. The top 30 genes from the highest module (magenta) were visualized with VisANT 5.0 (Hu et al., 2008) for gene connections with the thresholds of 0.02 and the result figure was shown in Figure 4D.

The Ingenuity's upstream regulator analysis in the Ingenuity Pathway Analysis (IPA, Ingenuity Systems, and www.ingenuity.com) was used to identify the upstream transcriptional regulators. Upstream regulators are identified based on knowledge of expected effects between transcriptional regulators and their target genes and, thus, can explain the observed expression changes in our gene dataset. During the IPA analysis, analysis settings included a cutoff for expression p-value (set at 0.01), reference set of Ingenuity Knowledge Base (Gene Only) and inclusion of both direct and indirect relationships. The analysis used Fisher's Exact Test to calculate p-values.

Alternative splicing events (ASE) were identified by the software Splicing Express (Kroll et al., 2015). Briefly, All GTF files created by Cufflinks (Roberts et al., 2011) were inserted with the mouse genome reference sequence (mm10) to run the software Splicing Express that is able to identify patterns of simple ASE, including exon skipping, alternative 5′ and 3′ splicing borders and intron retention (Roberts et al., 2011). Because of the software limitation, ASE results from only 122 genes that were selected and identified from the 15 most enriched GO biological process terms were statistically analyzed.

Student t test was used to detect a statistical significance in MWM performance between the high and low MWM performance. Two way repeated measures ANOVA (group × MWM days) was conducted to test MWM latency in each day followed by the post-hoc test of Holm-Sidak method. Normalized gene expression data from DESeq were analyzed by DESeq, PANTHER database, and WGCNA packages, as described above. Briefly, Differential expression (DE) data were generated using DESeq, to investigate significant differences between the groups. Corrected levels of significance and q-values were calculated to control the false discovery rate for multiple tests based on the Benjamini-Hochberg equations (Benjamini and Hochberg, 1995).

Normalized gene expression data from DESeq were statistically analyzed by WGCNA that computed all correlation and significance. Pearson product-moment correlation coefficient r was calculated between normalized gene expression data and individual rodents' MWM performance. Hierarchical cluster analysis with a Ward method was used to cluster genes and to provide a summary of selected gene sets. Gene lists from the DE analysis and correlation analysis were analyzed with PANTHER database. In addition, PANTHER statistical over- and under-representation test, by which expression levels of genes in each cluster are statistically tested to be correlated with categories of biological process or molecular function, was performed with gene lists and Log₂FC data. Best subset regression analysis was used with a selected gene set to identify best regression models. Bland-Altman plot (Bland and Altman, 1986) analysis was used to compare predictive values from regression models with true measured values from RNA-sequencing data and, thus, to evaluate an agreement between the regression models and RNA-sequencing data. All analyses were conducted by using statistical software of JMP® Pro (ver. 11.2.1, SAS Institute Inc.), SigmaPlot (ver. 12.3, Systat Software Inc.) and Microsoft Excel 2013. Alternative splicing event data from the Splicing Express (Kroll et al., 2015) were analyzed by Chi-squared and Fisher's exact tests in R (version 3.2.4) to test whether differences of the occurrence in ASE types between the high and low MWM performance groups were significant. Data are represented as mean ± standard error of mean (SEM), and an overall outline for data analyses is shown in Supplementary Figure S1.

Twenty-eight mice completed MWM testing over 5 days. The average time to locate the platform across the 5 days was determined (Figure 2A; MWM performance: mean = 24.67, median = 23.82, SD = 7.62, range: 26.55, CV = 30.90), and utilized to group animals based upon the distribution of their performance. As shown in Figure 2B, there was a significant difference between the means of the high performing mice (n = 9) and the low performing mice (n = 9). MWM latency in each day also resulted in significant group difference (Figure 2C), but no interaction between group and MWM days was detected. There was a significant group difference in the average of total distance moved during the 5 MWM days (Figure 2D), but no significant difference in the distance moved during the probe trial was detected (Figure 2E). The MWM latency was significantly correlated with total distance moved during the 5 MWM days (r² = 0.70, p < 0.0001), but no correlations of total distance during probe trial with MWM latency and total distance during the 5 MWM days were detected by Pearson product moment correlation analysis. No statistical group differences were observed on body weight at the time point of tissue collection (5 weeks), average body weight (1–5 weeks) and body weight change from the 1st week to the last week of the experiment (Supplementary Table S1).

RNA from hippocampal samples from all mice were sequenced. Quantified gene expression data (37,991 genes) from HTSeq analysis were normalized by DESeq, and, then, filtered based on coverage. The distribution of fold changes with p-values for these normalized-and-filtered genes (16, 565 genes) is depicted in the volcano plot (Figure 2F). Among them, 521 genes (red-colored dots) resulted in an absolute r (|r|) value greater or equal to 0.5 computed between gene expression and MWM performance. Most of these genes with |r| ≥ 0.5 were implicated in significant group differences in MWM performance based on a typical threshold of significance for RNA-sequencing data (Love et al., 2014).

The list of significant DE genes (n = 1,250) were entered into the website of PANTHER pathways analysis, and 111 pathways were identified (Figure 3 and Supplementary Table S2). The PANTHER statistical over- and under-representation test was also performed. The low MWM performance group showed under-expression in some intracellular and neurotransmitter receptor signaling pathways, including metabotropic glutamate, nicotinic acetylcholine, dopamine, corticotropin releasing factor receptor, GABA-B receptor II signaling pathways. The pathways overexpressed by the low MWM performance group included inflammation mediated by chemokine and cytokine signaling pathway, the gonadotropin releasing hormone receptor pathway, wnt signaling pathway, integrin signaling pathway, PDGF signaling pathway, and Alzheimer's disease. The most enriched GO biological processes (Tables 1, 2) that were detected by the PANTHER statistical over- and under-representation test include different receptor signaling pathways, apoptotic-related signaling pathways, immune related pathways and metabolism related pathways. Furthermore, the low MWM performance group showed statistically significant over expression in negative regulation of cell projection organization, negative regulation of neuron projection development, and innate immune response. The statistically significant under represented GO biological processes by the low MWM performance group included, but were not limited to, regulation of molecular function, regulation of apoptotic process, and negative regulation of nucleic acid-templated transcription.

To evaluate the differential expression of RNA profiles with MWM performance, we conducted weighted gene-coexpression correlation network analyses (WGCNA). Twelve modules (groups of genes) were merged into 10 based upon the similarity of their expression profiles, known as eigengenes (Figures 4A,B). When the associations between individual genes within each module and the trait (MWM performance) used in this study were quantified, a statistical significance was detected between the modules and the trait of rodent's MWM performance (Figure 4B). All 10 modules resulted in significant correlations with the individual MWM performance: we found positive and inverse correlations with 6 and 4 modules, respectively. The 3 highest correlations were detected by the module-trait association analysis for the magenta (r = 0.74, p = 0.0004), pink (r = 0.67, p = 0.002) and turquoise (r = −0.67, p = 0.003) modules. The meta-modules association analysis, one of WGCNA functions, shows the strength of mutual correlations between modules and the trait, identifying which module has the strongest mutual correlation with the trait. The magenta module showed the strongest mutual correlation with individual MWM performance (Figure 4C). The top 30 connected genes in the magenta module were depicted (Figure 4D), and the plot shows network connections whose topological overlap is above the thresholds of 0.02.

GO enrichment analysis was performed and the 20 best GO terms for each of 10 modules were identified (all listed on Supplementary Table S3). Briefly, among the significant GO terms, the biological process “innate immune response” was detected for both greenyellow and magenta modules. Four different receptor signaling pathways were detected from 3 different modules. GO terms relative to “splice” were detected for the brown and gray modules. The most significant GO terms include granulocyte differentiation (GO:0030851, enrichment p < 0.001, # of gene in magenta module = 3), protein targeting to peroxisome (GO:0006625, enrichment p < 0.001, # of gene in pink module = 3), mitochondrial transport (GO:0006839, enrichment p < 0.001, # of gene in turquoise module = 8), regulation of nitrogen compound metabolic process (GO:0051171, enrichment p < 0.001, # of gene in brown module = 32), neuron part (GO:0097458, enrichment p = 0.0018, # of gene in turquoise module = 29), and inflammatory response (GO:0006954, enrichment p = 0.002, # of gene in turquoise module = 10). The largest number of gene in a module was detected from the blue module: biological process (GO:0008150, enrichment p = 0.002, # of gene in turquoise module = 329).

Individuals' MWM performance during 5 days was significantly correlated with the expression of 411 genes (p < 0.050 and |r| ≥ 0.5, Table 3). The PANTHER statistical over- and under-representation test was performed, and many metabolic processes were found to be significantly over expressed in the low MWM performance group (Table 4). Interestingly, 5 processes detected are related to RNA regulation: mRNA processing (GO:0006397), RNA metabolic process (GO:0016070), RNA splicing via transesterification reactions (GO:0000375), mRNA splicing via spliceosome (GO:0000398) and RNA splicing (GO:0008380). These 5 RNA regulation-related processes are all significantly overexpressed in the low MWM performance group. Figure 5 shows biological processes with the 1st level of sub-categories, and the top biological processes are similar to the ones from previous analyses: such as, metabolism, cell communication, biological regulation, transport system, and immune system. The hierarchical cluster analysis with the highly correlated genes (p ≤ 0.010) resulted in a clear group separation and significant differences in the expression levels of these 82 genes (Supplementary Figure S2).

To identify potential transcriptional regulators that can explain the expression changes in our genes between the low and high MWM performance groups, the Ingenuity's upstream regulator analysis was conducted. The Ingenuity's upstream regulator analysis identified 16 different molecule types and 247 upstream transcriptional regulators (all listed on Supplementary Table S4). Among them, 7 cytokines (IL1B, IL3, IL5, LIF, MIF, TNF, and TNFSF13B) were identified and the reduced activation was observed from IL1B (activation z score = −1.954, p for overlap = 0.041). Nicotinic acetylcholine receptor was also identified and negative activation z score (−1.997) was calculated for this (p for overlap = 0.001). Thirty four transcriptional regulators were identified, but activation z scores for only CREB1 (−1.969) and CREM (−1.982) were identified with significant p-values for overlap (0.002 and 0.001, respectively). Network for upstream regulators was created (Figure 6) and 28 regulators were identified to be associated with learning behavior (outlined with pink color). Among these 28 regulators, CREB1 and CREM regulators were predicted as inhibition and TP63 and ICOSLG regulators were predicted as activation. When overlaid with biological functions, 15 upstream regulators were significantly associated with learning behavior (p = 2.43E-11), and 23 upstream regulators were significantly associated with nervous system development and function (p-value range: 1.32E-7–6.51E-7). When overlaid with diseases, neurological diseases, including Alzheimer disease (p = 5.02E-7), were identified to be significantly associated with 38 upstream regulators (p-value range: 1.4E-11–5.13E-6). Additionally, the network showed that FOS, NR4A2, WNT6, and DUSP1 regulators were highly decreased in our gene dataset. With our gene dataset, the IPA identified 5 networks (score ≥ 10, Table 5). The network 1 contained 17 upstream transcriptional regulators (Figure 6). Some upstream transcriptional regulators were also identified from the networks 2 (NCL), 3 (TRIB1). 4 (DIAPH3, GRHL2, MCIDAS, TP63), and 5 (CAPN10, CHRNA5, CTHRC1, RYR2).

Because unique biological processes identified a potential contribution of RNA splicing to group differences (Table 4), analyses of ASE were carried out using the recently developed software Splicing Express (Kroll et al., 2015). No statistical differences were detected in the number of Genes (X-squared = 0.81714, p = 0.85) and the number of Total Events (X-squared = 1.4336, p = 0.70) between the high and low MWM performance groups (Supplementary Table S5). Because of the difficulty in pulling out all ASE data from the Splicing Express, 122 genes in the 15 most enriched GO biological process terms (Supplementary Table S6) were chosen to analyze ASE. Overall, similar patterns in the ASE types between the groups were discovered (Supplementary Figure S3). However, we found group differences in the ASE exon skipping type in genes (Figure 7A, p-value from Fisher's exact test was 0.043). The low MWM performance group has more genes with 4, 6 and 10 exon skipping events but fewer genes with 5, 7, and 8 exon skipping events. Alternative 3′ and 5′splice sites and intron retention for these genes are depicted on Figures 7B–D (p-values from Fisher's exact test were 0.82, 0.32, and 0.12, respectively). Because the Figure 7D showed a visually greater difference in certain intron retention events (especially 2, 3, 4, and 5 events), the Fisher's exact test was run with only 1 through 5 intron retention events, then a statistical significance was detected (p = 0.048). These ASE data analyses suggest that the number and types of ASE for some genes were significantly different between the low and high MWM performance groups.

To analyze whether individual's MWW performance could be predicted from the expression of selected genes, best subset regression was performed with the top highly correlated genes. This allows the identification of the best regression model for cognitive performance under chronic stress. From the top 16 genes, the best subset regression method generated 10 potential regression models (from a total of 16 models, Supplementary Table S7). Variance inflation factors (VIFs) of each model were computed and excluded the models 5–10 due to VIFs > 2 although these models presented higher r²-values than the remaining 4 models. The 4 regression models (0.68 ≤ r² ≤ 0.96) chosen are depicted with their statistical information in Figures 8A,B. To validate these models, we added additional whole transcriptome RNA-sequencing data from 7 different strains (average performing strains) that were excluded from grouping, and drew a Bland-Altman (BA) Plot to analyze the agreement between predicted values (from the equations obtained from the 4 regression models) and measured values (from normalized DESeq data). As shown on Figures 8C,D, bias ranged from 1.18 to 2.41, and at least one sample from all models was detected above the 95% upper limits of agreement CI. Among these models, the model# 4 showed the smallest bias (bias = 1.19, 95% Bias CI: −0.70 to 3.08) with the smallest range of limits of agreement (−7.76, 10.14), and the differences of 24 samples (out of 25 samples; thus, 96%) between the prediction and measurement values were positioned between the 95% upper and lower limits of agreements.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.