8-week-old male C57BL/6J mice purchased from The Jackson Laboratory (Bar Harbor, ME, United States) were used for all experiments. Upon arrival, mice were group housed, 3-4 mice per cage, under a 12-h light/12-h dark cycle with ad libitum access to food and water. Testing and handling was performed during the light phase of the cycle. Behavioral tests and network analyses were conducted using groups of n = 10. Mice from the Morris Water Maze training group and the cage control group were co-housed, and all mice received equal handling throughout the study. All procedures were conducted in accordance with protocols approved by the University of Calgary, Health Sciences Animal Care Committee, following the guidelines of the Canadian Council for Animal Care.

In order to provide mice with chronic cognitive stimulation, we trained half of the mice on a repeated acquisition and performance testing variant of the Morris Water Maze (MWM; see Figure 1A for an outline of the testing schedule and Figure 1B for the maze itself) (Spanswick et al., 2007). This task was chosen because it provides considerable flexibility in design, whereby a hidden platform can be moved to many different locations within the maze to encourage continuous learning. Furthermore, the spatial learning which occurs during MWM training is supported by many of the same neuroanatomical regions as contextual conditioning (Jo et al., 2007; Miller et al., 2014; Giustino and Maren, 2015; Kwapis et al., 2015; Milczarek et al., 2018). In this version of the task, the hidden escape platform is moved every second day which requires mice to repeatedly acquire new spatial memories throughout a 10-week training period. We used 10 different platform locations, with each platform location occurring on 2 separate occasions. Mice were trained 4 days per week (2 platform locations). During each daily session, mice were given four trials. Each trial lasted a maximum of 60 s and was initiated by placing the mouse gently into the pool, facing the wall. The start location was from a different cardinal compass position around the pool for each trial and the order of start locations was randomized each day. Trials were terminated once the mouse located the hidden platform. If the platform was not found after 60 s, mice were gently guided to the platform by the experimenter. Once on the platform, mice were given 15 s to remain on the platform before being returned to their cage. Trials were interleaved, whereby each mouse performed their first trial before the first mouse performed its second of four daily trials. This resulted in an intertrial interval of approximately 10 min. The circular pool had a diameter of 150 cm and a depth of 50 cm. The pool was filled so that the water level was 2 cm above the surface of a circular escape platform that had a diameter of 11 cm. The water was made opaque using white non-toxic tempera paint. The water was kept at a constant temperature of 22°C and stirred and cleaned of debris before each trial. Automated tracking software (ANY-Maze, Stoelting, Wood Dale, IL, United States) was used to record and analyze swim behaviors in the pool, primarily the distance traveled prior to locating the hidden platform. When analyzing these results, linear regression was applied to the mean distance traveled by each mouse across each training session to assess the influence of the memory of previous platform locations on the ability to learn new platform locations. The extent to which mice learned platform position within blocks of consistent locations was assessed by examining the mean slope of the regression lines between distance traveled at all first-day trials and between all second day trials.

After the conclusion of the Morris Water Maze training protocol, all mice were trained in contextual fear conditioning. Training was conducted in sound-attenuated chambers with grated floors through which shocks (0.5 mA; 2 s) were delivered (Ugo Basile, Gemonio, Italy). Mice were first allowed to acclimate to the chamber for 2 min prior to the presentation of a series of 3 shocks, each separated by an interval of 90 s. 24 h after the training session, mice were returned to the conditioning chambers for a 6-min retention test. During this test, no shocks were administered, and behavior was monitored via an overhead infrared camera in conjunction with an automated tracking software (ANY-Maze, Stoelting, Wood Dale, IL, United States). The chamber was cleaned using 70% ethanol and allowed to dry before and after each trial. When analyzing these results, freezing criteria was defined as bouts of a minimum of two seconds without ambulation. The percentage of the trial spent exhibiting freezing behavior was compared between groups.

Mice were transcardially perfused with 0.1 M phosphate buffered saline (PBS) followed by 4% formaldehyde 90 min after retention testing. Brains were then extracted and post-fixed in 4% formaldehyde for 24 h. Fixed brains were cryoprotected in 30% W/V sucrose solution at 4°C until no longer buoyant. From cryoprotected brains, serial coronal sections with a thickness of 40 μm were cut on a cryostat (Leica Biosystems, Concord, ON, Canada) and stored in 12 series at −20°C in antifreeze solution.

When conducting labeling for c-Fos expression, all tissue was processed concurrently. Tissue sections were washed 3 times (10 min per wash) in 0.1M PBS before being incubated in a primary antibody solution of 1:2000 rabbit anti-c-Fos primary antibody (226 003, Synaptic Systems, Göttingen, Germany), 3% normal donkey serum, and 0.03% Triton-X100 for 48 h at room temperature on a tissue shaker. Tissue sections were washed 3 × 10 min in 0.1M PBS before secondary antibody incubation. The secondary antibody solution was composed of 1:500 donkey anti-rabbit Alexa Fluor 488 (111-545-003, Cedar Lane Labs, Burlington, ON, Canada) in PBS for 24 h at room temperature. Sections were then transferred to 1:2000 DAPI solution for 15 min before being washed 3 × 10 minutes in 0.1M PBS. Labeled sections were mounted to glass slides and coverslipped with PVA-DABCO mounting medium.

Quantification of fluorescent c-Fos labeled cells was conducted using a custom semi-automated segmentation and registration pipeline (Figure 2A). All slides were imaged as a single batch using an Olympus VS120-L100-W slide scanning microscope (Richmond Hill, ON, Canada). Images were collected using a 10x objective with a numerical aperture of 0.40 and a Hamamatsu ORCA-Flash4.0 camera. Labeled c-Fos was imaged using a FITC filter cube and a 9.00 V lamp at an intensity of 100% and an exposure time of 140 ms. DAPI staining was imaged under the same conditions, but with a DAPI filter cube and an exposure time of 65 ms. Cells expressing a c-Fos label were segmented using the machine learning-based pixel and object classification program, Ilastik (Berg et al., 2019). To further prepare Ilastik output images and DAPI channel photomicrographs for regional registration, a custom plug-in was written for ImageJ. The pixel intensity threshold of the Ilastik outputs was adjusted so as to only contain objects which the program determined to be within the correct range of pixel intensities and shapes. To compensate for inadequate regional area measurements at an image-by-image level in the atlas registration software, a mask of evenly spaced binary points was generated from the DAPI channel image. The pixel intensity thresholds of these images were adjusted to create a binary mask in the shape of the tissue section. Grid lines were then overlaid to create a mask of binary points arranged in a square grid in the shape of the tissue section. Adjacent binary points were spaced by 22 μm, therefore, each point in the mask accounted for an area of 484 μm².

Next, tissue sections were registered to plates of the Allen Mouse Brain Atlas using the R-based Whole Brain software (Fürth et al., 2018). Using this software, DAPI channel images were used as references to which the atlas plates were aligned. The number of segmented c-Fos labeled cells per neuroanatomical region was quantified in Whole Brain. Similarly, the binary point masks were processed to count the number of points in each region. Regional areas were then approximated using a Cavalieri-based point counting approach, whereby the number of mask points in each region was multiplied by the area accounted for by each point. This allowed for the c-Fos labeled cells to be normalized by area and presented as regional cell densities.

A separate cohort of mice was used for the validation of c-Fos labeled cell segmentation and regional area approximation. c-Fos immunostaining and imaging was identical to the methods described previously. To generate ground truth cell counts as a gold standard for our automated counting procedure, ten 500 μm x 500 μm regions of interest (ROI) were randomly generated from each of several regions including the basolateral amygdala, CA1, dentate gyrus, paraventricular nucleus, and the retrosplenial cortex. These ROIs were processed through our Ilastik label segmentation pipeline (representative raw and processed images in Figure 2B). In addition, the same ROIs were hand counted independently by 4 experimenters blind to the automated cell count results. The total numbers of cells counted across all counting boxes were then compared to assess whether or not Ilastik could segment fluorescent c-Fos labels within natural and acceptable inter-rater variability.

To assess the accuracy of regional area approximations, areas generated using the Cavalieri-based point counting approach were compared to the areas of these same regions which were manually traced in ImageJ. During this analysis, five photomicrographs of each of the following regions were examined: basolateral amygdala, CA1, dentate gyrus, paraventricular nucleus, and the retrosplenial cortex.

We focused on a selection of 97 regions based on our ability to discriminate these regions using a DAPI stained image as reference (see Supplementary Table 1 for list of regions and abbreviations). From this list of regions, c-Fos label densities were cross correlated within each group to generate pairwise correlation matrices. Correlations were filtered by statistical significance and a false discovery rate of 5% (Benjamini and Hochberg, 1995; Bassett et al., 2009). For network analyses, correlation matrices were binarized to adjacency matrices based on Pearson’s correlation coefficient and statistical significance (r > 0.9; α = 0.005). This threshold allowed for sufficient network density to study global brain dynamics, while still limiting the analyses to only the strongest and most biologically plausible connections (Schneidman et al., 2006). To ensure that the thresholding parameters did not bias network analyses, additional adjacency matrices were generated using either more (r > 0.95; α = 0.0005) or less (r > 0.8; α = 0.05) conservative thresholds. To analyze adjacency matrices as network graphs, the 97 neuroanatomical regions were plotted as nodes. Connections were drawn between nodes whereby correlations surpassed correlation matrix thresholding parameters.

Graph theoretical analyses were applied to network graphs to examine global and local properties of the network. These analyses were guided by the use of the Brain Connectivity Toolkit (Rubinov and Sporns, 2010), the SBEToolbox (Konganti et al., 2013), and other custom analyses. Network properties examined include node degree, network density, global efficiency, betweenness centrality, Katz centrality, and network resiliency. In the following definitions, N is the array of nodes in the network represented by adjacency matrix A. The number of nodes in the network is represented by n and the number of connections between nodes is l. The variable a_ij is the index into the adjacency matrix which indicates the connection status of nodes i and j. The presence of a connection is represented by a_ij≠ 0 (Rubinov and Sporns, 2010). Node degree is the number of connections that link a node to the rest of the network (Rubinov and Sporns, 2010).

Network density is a metric of network dispersion. It is expressed as a proportion of the number of connections in a given network over the number of connections which would be required to saturate a network of the same size (Rubinov and Sporns, 2010).

Global efficiency is defined as the inverse of the average shortest path of connections between all possible pairs of nodes (Latora and Marchiori, 2001; Achard and Bullmore, 2007; Rubinov and Sporns, 2010).

Betweenness centrality and Katz centrality are measures which can be used to assess the importance of a node in the effective communication of a network. Betweenness centrality quantifies the number of shortest paths between nodes that pass-through a given node (Freeman, 1978; Brandes, 2001; Rubinov and Sporns, 2010).

Katz centrality applies an eigenvector approach to this metric by weighting the connections involving more highly connected nodes more heavily than those from lesser connected nodes when considering the makeup of the shortest paths which pass through a given node (Katz, 1953; Hubbell, 1965). The attenuation factor, α, used for this analysis was 0.1 (Zhan et al., 2017).

Network resiliency was assessed through targeted node deletion and an assessment of the size of the largest community of connected nodes and global network efficiency with each deletion. Nodes were targeted for deletion in decreasing order, from nodes with the highest degree to those with the lowest. Degree was recalculated after each deletion and the list was reordered accordingly.

Network metrics were both compared across conditions as well as used to assess small world-like network properties compared to random control network topology. Small world network distribution can be described as being efficient at both local and global scales (Watts and Strogatz, 1998). Random null control networks were generated for both the cognitive training and control groups and were matched for network size, overall degree, and degree distribution. Local efficiency was assessed by comparing mean clustering coefficients, while global efficiency was assessed by comparing bootstrapped global efficiency values with one hundred replacements. Networks were considered to display small world-like properties if they had displayed both the high global efficiency characteristic of a random network and increased local efficiency relative to random networks (Wheeler et al., 2013).

Behavioral data from all tasks was recorded and analyzed using ANY-Maze (Stoelting, Wood Dale, IL, United States). All t-tests and Two-Way ANOVA for comparing behavioral data, regional c-Fos expression, segmentation and regional area approximation validation, and functional connectivity networks were conducted using GraphPad Prism (GraphPad Software, San Diego, CA, United States). GraphPad Prism was also used to conduct the linear regression used to assess inter-position memory and repeated new learning in the Morris Water Task. The analysis of functional connectivity networks was conducted using MATLAB. Figures were generated using MATLAB, Cytoscape, and GraphPad Prism.

A repeated acquisition and performance testing variant of the Morris Water Maze was used to provide chronic cognitive stimulation (Figures 1A–C). Simple linear regression was applied to assess inter-position memory and the presence of repeated new learning. The slope of the best fit line of the linear regression applied to the mean distance traveled by each mouse across each training session was −0.067. When examining the mean distances traveled across all first days at a given platform position, the line of best fit had a slope of −0.068. Across all second days at a given platform position, this slope was determined to be −0.069. Within each platform location, the line of best fit of the linear regression yielded a mean slope of −0.97 which differed significantly from zero and was indicative of improved performance over time within the same platform location (one sample t test; p = 0.0010).

To assess the generalization of improved cognitive performance following long-term spatial learning, mice were trained and tested in a contextual fear conditioning paradigm (Figure 1D). The percentage of time that mice exhibited freezing behavior was compared between groups and across the training and retention test sessions. Across both sessions, there were no significant differences in freezing behavior (Figure 1E). However, when the retention was divided into a first half and a second half, mice who had received cognitive training displayed increased freezing behavior during the first half of the test and then decreased freezing during the second half of the test (Figure 1F).

To assess the reliability of the semi-automated c-Fos segmentation and mouse brain atlas registration pipeline used in this study (Figure 2A), we compared c-Fos counts obtained using this pipeline to those gathered manually (see representative images Figure 2B). We found that the number of c-Fos labeled cells quantified using Ilastik processing fell within the range of values counted manually by four different experimenters across a subset of regions with varying levels of background autofluorescence (Figure 2C). The inter-rater reliability was determined to be 15% across the datasets as a whole and the semi-automated cell counts were within 3.75% of the average of the manual cell counts. Manual and automated cell counts were highly correlated across the sampled brain regions (Figure 2D). WholeBrain region registration also produced regional area approximations that correlated highly with hand-traced regional area values (Figure 2E).

Changes in regional c-Fos expression density are depicted in Figure 2F as the percent change from the regional c-Fos expression density from the control condition to the group which had underwent long-term spatial learning. We observed decreased c-Fos expression density in all brain regions that were analyzed in mice that had received prior spatial training. With respect to brain-wide activity, prior spatial training resulted in an overall significant decrease in c-Fos expression (Supplementary Figure S1. Unpaired t test; p = 0.0003).

Analyses of cross-correlated regional c-Fos expression density revealed differences in global functional connectivity network topology. On a global scale, we observed a reorganization of connections throughout the brain (Figures 3A–F). Relative to control conditions, mice that underwent prior cognitive training exhibited an increase in the overall density of functional connections during subsequent contextual fear memory retrieval (Figure 3G). Furthermore, the organization of these networks after prior cognitive training resulted in increased global efficiency relative to the control condition (Figure 3H). These increases were also present in the networks constructed with both more or less conservative thresholds, indicating that this effect was not an artifact of the thresholding level (Supplementary Figure S2).

Networks from both control and spatial learning groups exhibited heavy-tailed degree distribution characteristic of a small-world network, with the majority of nodes making very few connections and a lesser number of nodes carrying disproportionate importance to the overall connectivity of the network (Figure 4A) (Bassett et al., 2006; Bullmore and Sporns, 2009). Comparisons to random null networks also highlighted that both the control and spatial learning networks maintained the high global efficiency characteristic of random networks (Figure 4B) while displaying increased clustering (Figure 4C). Together, these analyses indicate that the functional connectivity networks engaged during memory recall in both the control and spatial learning conditions exhibit properties that are consistent with small-world topology.

Changes in network topology were also observed at the local level. The organization of local communities within global networks changes with long-term spatial learning. While the size of the giant component (GC) (Figures 5A–C) underwent very little change with MWM training, differences arose in the connectivity patterns within this component. Within the GC, MWM training increases the mean number of connections per node (Figure 5D). The changes in connectivity coincided with changes in network resiliency. When faced targeted deletion of nodes, with deletions occurring in the order of decreasing degree, long-term spatial learning increased the ability of the network to preserve its giant component size (Figure 5E) and global efficiency (Figure 5F).

Coinciding with these changes in network resiliency were changes in Katz centrality. Katz centrality is a measure of centrality which differentially weighs connections based on the degrees of the nodes involved (Katz, 1953). This measure has previously been shown to correlate highly with neuronal activity compared to other measures of centrality (Fletcher and Wennekers, 2018). While most nodes in the control network had similar Katz centrality vectors (Figure 5G), MWM training increased the centrality of a subset of regions (Figure 5H). There was considerable overlap between these regions with increased Katz centrality and the regions in the most densely connected region of the GC. This was further corroborated by analysis of regional degree distribution (Supplementary Figure S3) and change in Katz centrality (Supplementary Figure S4) which further highlighted an increase in connectivity and of importance of numerous amygdala subregions following spatial learning.