Male and female C57BL/6 J (B6) mice aged 10, 30, or 60 weeks old (wo) at the time of arrival were allowed a week of acclimation to the vivarium prior to the start of experiments (N = 30, total; n = 10 per age group, 50% females; The Jackson Laboratory, Bar Harbor, ME). Female mice initially weighed (in grams; mean ± standard deviation) 19 ± 0.93, 32 ± 4.6, 29 ± 2.8 and males 26.7 ± 2, 30 ± 2.5, 46.8 ± 5.6 (for 10, 30, 60wo, respectively). The selected age groups represent 3 stages of adulthood. B6 mice at 3–6MO (12–24wo) are considered mature adults and at 10–14MO (40–56wo) are considered middle age (Flurkey et al., 2007). The youngest 10wo group (2.5MO) was tested just before mature adulthood and the 60wo group was just at the end of middle age (15MO). At 30wo (7.5MO), mice start to show early signs of senescence (such as retirement from breeding at 8MO) and were included here as reference to compare early adulthood and late middle age changes (Flurkey et al., 2007). Mice were housed in age- and sex-matched groups of 3–4 in a temperature- and humidity-controlled room, inside conventional air filtered cages (dimensions: 29 × 18 × 13 cm) with food and water available ad libitum (vivarium lights on from 07:00 to 19:00 h). All procedures received prior approval from the Institutional Animal Care and Use Committee of the University of Florida and follow all applicable NIH guidelines.

Figure 1 shows the experimental timeline. All behavioral tests and functional neuroimaging were conducted over an 8-week period. Behavioral tests included locomotor activity in a novel environment, sociability and social recognition tests, sucrose preference tests, and cued fear conditioning (Figure 1). Functional magnetic resonance imaging (fMRI) was conducted a week after sucrose preference testing and 2 weeks before the final behavioral test (fear conditioning). At the end of the experiments, mice were deeply anesthetized with 5% isoflurane, and brain tissue harvested for the metabolomic assay of excised prefrontal cortex and hippocampus tissue. Imaging was thus carried out when mice were 14, 34 and 63wo and these were euthanized 8 weeks after arrival, at 18, 38 and 68wo. However, to simplify presentation of the results, we keep group names in figures as 10, 30 and 60wo.

An automated rodent activity monitoring cage system (Accuscan Instruments, Columbus, OH) was used to measure spontaneous locomotor activity in a novel environment (Bruijnzeel et al., 2016). The day before locomotor measurements, mice were acclimatized for 1 h to the testing room. The activity monitors are made of clear Plexiglas (40 cm³) with an outer steel frame containing infrared emitting and receiving panels. Each panel has 16 equally spaced (2.5 cm) infrared beams crossing the length and width of the cages at a height of 2 cm (horizontal beams). Several patterns of sequential and repetitive infrared beam interruptions were registered and sent to a programmable counter, which sorted beam counts and sent sorted data to a PC running Versamax^™ software. The sorted activity counts were classified as sequential beam breaks (forward locomotion), repetitive beam breaks (putative repetitive behaviors), total, center and margin distances traveled (in cm). Activity was measured during a single 30-min test.

The sucrose preference test was conducted as described previously, with modifications (Bruijnzeel et al., 2016; Chellian et al., 2020). A decrease in preference to consume a 2% sucrose solution over water suggests anhedonia (D’Aquila et al., 1997). ‘Leakproof’ bottles (100 mL, graduated by 1 mL) with stoppers and sipper tubes were purchased (Braintree Scientific Inc., Braintree, MA; model no. WTRBTL). Filtered home cage lids with two bottle openings (Allentown, Inc., Allentown, NJ), were used for two bottle fluid delivery in mouse home cages. Mice were temporarily single housed for testing and once testing concluded they were returned to their original cages. First, mice were acclimatized to bottles filled with 2% sucrose solution for 3 days (96 h). On each of the 3 days, bottles were weighed (in grams) to monitor fluid intake and replaced with a full bottle. The 96 h acclimation period was followed by a single 24 h two-bottle choice test in which one bottle was filled with 2% sucrose solution and the other with water. The side of the cage lid where the sucrose and water filled bottles were placed was counterbalanced across animals. Water and sucrose consumed was measured after a 24 h period (McCall et al., 2013). Bottles were weighted and values converted from grams to ml. An index (S_I) of how much sucrose (s) consumed vs. water (w) was calculated as: si=s−ws+w.

We used a modified version of a previously reported sociability and social recognition test protocol (Tan et al., 2019). Test arenas were made of clear Plexiglas, divided by walls into 3 equal sized quadrants (arena total length x width x height dimensions: 60 × 45 × 20 cm). The middle quadrant is the neutral zone, and the side quadrants (left and right) were used for the placement of ‘bait’ mice for sociability and recognition tests. Each side quadrant can be entered through a central door in the dividing walls. A pair of clear plastic ‘bait’ cups covered with equally sized and equidistant perforations for air flow were inverted over the ‘bait’ mouse for sociability testing. Test mice were acclimatized to the room and the neutral quadrant of the test arena on 3 days, 15 min a day. During sociability testing, a novel (unfamiliar) mouse (an adult mouse of the same strain and sex; age range 3–6mo) was placed under a cup and the entrance to the ‘bait’ side of the arena was opened for exploration by the test mouse while the other side is kept sealed. In preliminary tests, this encouraged mice to explore the unfamiliar ‘bait’ mouse more vigorously during a first encounter. Mice were tested for 10 min for time spent exploring the unfamiliar mouse. Upon completion of the first test, and after a 5-min break to clean the arena, mice were placed back into the arena in the presence of a novel mouse and the previously encountered (now familiar) mouse and re-tested for 10 more minutes with both sides opened for exploration. In the first test, the degree of exploration of an unfamiliar mouse (“novel 1”) is considered here as a measure of the degree of sociability. In the second test, less exploration of the familiar mouse (“familiar 1”) relative to the new unfamiliar mouse (“novel 2”) is a sign of social recognition (Crawley et al., 2007). Behaviors were recorded using a digital camera directly connected via FireWire to a Linux computer. We used the ‘dvgrab’ tool from the Unix command terminal to program video capture settings and store the recorded videos.

Behavioral analyses were carried out using DeepLabCut (DLC, version 2.2.1; Mathis et al., 2018; Nath et al., 2019). We labeled 120 frames of mouse trajectories within the test cages taken from 6 different videos and 95% of those labels were used to train a machine learning algorithm. A ResNet-50-based neural network with default parameters for 2 training iterations was used for training. Data (digital videos) were split into training and test datasets. The training dataset was used to train a DLC model. Conversely, the test dataset was used to evaluate the performance of the model. We validated the DLC model with 7 train-test-split shuffles. The model with the smallest mean average Euclidean error (MAE) was selected. Before model selection, we applied a p-value cutoff of 0.6 to exclude model predictions with a likelihood under 0.6. The final selected model had a MAE of 2.9 pixels. This neural network was used to analyze videos from similar experimental settings to generate trajectory data. Using the R packages ‘trajr’ (McLean and Skowron Volponi, 2018) and ‘DLCanalyzer’ (Sturman et al., 2020), spurious trajectories were removed before analyses. The videos were analyzed in two parts, the first sociability test (T1) and the second social recognition test (T2). The immediate area surrounding and including the mouse ‘bait’ cup on one side of the arena was designated ‘ROI A’ and the other ‘ROI B’ (region of interest A and B, respectively). A python script was compiled to generate the number of entrances, the velocity (in pixels per frame or ppf), and the duration of initial approach and time spent within the ROIs.

Fear conditioning was conducted 2 weeks following the fMRI studies. The fear conditioning cages were housed in sound attenuation chambers. The chambers had internal computer-controlled lighting, white noise generator, a fan, and control units for distinct aspects of the fear conditioning assay. These included audio speaker controls for tone generation, shock controllers for the internal cage floors, and infrared and white lighting controllers. The internal operant cage was made of translucent Plexiglas with a steel frame and steel shock grid floor connected to a computer controlled current stimulator (Med-Associates, Inc. St. Albans, VT). Activity within the fear conditioning cages was captured by a monochrome video camera and video recordings were optimized via externally controlled near infrared (NIR), white light and sound controllers. All cage accessories were controlled by PC running Video Freeze Software™. Camera lens brightness, gain, and shutter settings were verified and adjusted before collecting data and were kept constant in all mice. For each session, the detected motion-sensitive signals were calibrated with the mouse inside the cage and NIR lights on. The estimated motion index was set to the same threshold level (10 a.u.) in all mice prior to data exporting or analyses.

Videos were acquired at 30 frames per second (fps) and the minimum freeze duration was 30 frames. For training on the first day, a 180 s baseline was followed by four presentations of a 20 s tone (90 dB, 5 kHz, 50 ms risetime) with the presentation of moderate level current to the entire grid floor for 1 s (0.9 mA; inescapable shock) on the final second before the tone ended. This 90 dB tone is above the auditory brainstem response measured in C57BL/6 J mice, which declines with age and reaches 60–80 dB at 5–8 kHz (Henry, 2004; Yu et al., 2011; Youn et al., 2020). Each of the four tone-shock presentations were separated by 190 s and the entire session lasted 17 min. On days 2 and 3 (recall tests at 24 and 48 h of training day 1, respectively), the same protocol was carried out but in the absence of a shock. On day 3 (modified context test), the cage environment was modified by placing a plastic smooth white floor and diagonal walls, which covered the grid floor, light, and speaker.

For each of the time epochs corresponding to the 4 tone-shock pairings, a 20 s pre-shock tone interval and a 40 s post-shock interval were used to quantify number of freezing bouts (freezing counts) in a total 60 s epoch (per day there were four 60 s epochs summed to provide total freezing counts per day). Motion-index was also used as a secondary measure of general locomotor activity inside the cage during the same 60 s epochs/day.

Images were collected on a magnetic resonance spectrometer tuned to 200.6 MHz proton resonant frequency (4.7 Tesla/33 cm actively shielded magnet). The MR system was equipped with Resonance Research Inc. spatial encoding gradients (RRI BFG-200/115-S14, maximum gradient strength of 670 mT/m at 300 Amps and a 120 μs risetime; RRI, Billerica, MA) and controlled by a VnmrJ 3.1 console (Agilent, Palo Alto, CA). A custom-built 8.8 cm x 11.6 cm radiofrequency (RF) transmit coil used in combination with a 3 cm x 3.5 cm actively decoupled/tunable RF receive quadrature surface coil (RF engineering lab, Advanced Magnetic Resonance Imaging and Spectroscopy Facility, Gainesville, FL).

Mice were scanned sedated under a continuous paranasal flow of 0.25% isoflurane gas (delivered at 0.04 L/min mixed with medical grade air containing 70% nitrogen and 30% oxygen, Airgas, Inc.) and a continuous subcutaneous infusion of dexmedetomidine. Prior to setup, mice were anesthetized under 2% isoflurane and administered an intraperitoneal (i.p.) injection of 0.1 mg/kg dexmedetomidine (at a 1 mL/kg volume). Isoflurane was reduced to 0.25% throughout the remaining imaging session. An infusion line was setup for subcutaneous delivery of dexmedetomidine over the course of scanning (0.1 mg/kg/ml at an infusion rate of 25 μL/h using a PHD-Ultra microinfusion pump, Harvard Apparatus, Holliston, MA). Functional MRI scans were collected at least 50 min after the i.p. injection. Respiratory rates were monitored continuously, and body temperature was maintained at 36–37°C using a warm water recirculation system (SA Instruments, Inc., New York).

We acquired one T2 weighted anatomical and an fMRI scan per mouse. The T2-weighted fast spin echo (FSE) sequence was acquired with the following parameters: effective echo time = 48 ms, repetition time (TR) = 2 s, echo train length = 8, number of averages = 12, field of view (FOV) of 19.2 mm x 19.2 mm and 0.8 mm thick slice, and a data matrix of 192 × 192 (0.1 μm² in plane) and 12 interleaved ascending coronal (axial) slices covering the entire brain from the rostral-most extent of the anterior prefrontal cortical surface, caudally toward the upper brainstem and cerebellum (acquisition time 9 min 36 s). Functional images were collected using a single-shot gradient echo planar imaging sequence with the following parameters: TE = 22 ms, flip angle = 90°, TR = 1.2 s, 600 repetitions, FOV = 19.2 × 19.2 mm and 0.8 mm thick slice, and a data matrix of 48 × 48 (0.4 μm² in plane) with 12 interleaved ascending coronal slices in the same position as the anatomical scan. Ten dummy EPI scans were run prior to acquiring data under steady state conditions. Respiratory rates, isoflurane and dexmedetomidine delivery, body temperature, lighting, and room conditions were kept constant across subjects.

The image processing workflow is illustrated in Figure 2 (Colon-Perez et al., 2019; Kotlarz et al., 2022; Sadaka et al., 2023; Sakthivel et al., 2023). Resting state processing was carried out using Analysis of Functional NeuroImages (AFNI; Cox, 1996), FMRIB Software Library (FSL; Jenkinson et al., 2002), and Advanced Normalization Tools (ANTs; Klein et al., 2009). Binary masks for anatomical and functional scans were created using ITKSNAP (Yushkevich et al., 2006). The brain binary masks were used for brain extraction prior to registration steps. Times series spikes were removed (3dDespike, AFNI), image repetition frames aligned to the first time series volume (3dvolreg, AFNI), and detrended (high pass filter <0.009 Hz using AFNI 3dTproject). Independent component analysis (ICA) using Multivariate Exploratory Optimized Decomposition into Independent Components (FSL MELODIC version 3.0) was used to assess structured ‘noise’ or artifact components in each subject scan, in their native space. Most, if not all ICA components in this first stage contained artifact signal voxels along brain edges, in ventricles, and large vessel regions. These components were suppressed using a soft (‘non-aggressive’) regression approach, as implemented in FMRIB Software Library (FSL 6.0.3) using fsl_regfilt (Jenkinson et al., 2002). A low-pass filter (>0.12 Hz) and spatial smoothing (0.4 mm FWHM) was next applied to the fMRI scans prior to registration steps. Post-regression ICA was carried out to verify removal of artifact components and preliminary assessment of putative ICA networks in individual scans.

Anatomical scans were cropped, and bias field corrected (N4BiasFieldCorrection, ANTs). Functional scans were cropped, and a temporal mean image was registered to the anatomical. The functional-to-anatomical registration matrix was applied to all fMRI volumes. A multi-T2 anatomical template was generated using ANTs ‘antsMultivariateTemplateConstruction2.sh’ script (cross correlation used as a similarity metric and a greedy SyN as transformation model). Individual subject-to-template transformation matrices were applied to fMRI datasets. A group ICA was applied to template-registered concatenated fMRI scans to generate a set of 60 nodes for ROI-based intrinsic fMRI signal extraction, pairwise cross-correlation statistics, and adjacency matrix construction. The ICA was set to an upper limit 60 component signal decomposition and a mixture model threshold of 0.7. The resulting 60 ICA components included peak signal (peak T statistic) nodes in midline prefrontal cortex (anterior cingulate), somatosensory regions, temporal, visual, auditory cortices, entorhinal cortex, several hippocampal nodes, subiculum, amygdala, thalamic nuclei, and striatum. Each of the 60 nodes were binarized and assessed to ensure non-overlapping coverage of brain areas on the template and registered scans (Figure 2). Both the template and ICA nodes were normalized to the Allen Brain Institutes common coordinate framework version 3 (CCFv3) atlas of the B6 mouse brain for final 3D visualizations (Lein et al., 2007). Signal timeseries were extracted from preprocessed fMRI scans with the assistance of 60 ROI mask overlays. This generated 60 individual ROI text files per subject that contained L2-normalized resting state signals as a vector of 600 data points.

The timeseries files were used in cross-correlations and in calculations of Pearson r coefficients for every pairwise combinations of ROIs (using corrcoef function in MATLAB). The resulting number of pairwise correlations was 1,710 per subject (after removing 60 self-correlations). Correlation coefficients were Fisher’s transform to ensure normality prior to statistical analyses. The center voxel coordinates for the ICA-based nodes normalized to the CCFv3 were used in 3D network visualizations using BrainNet viewer (Xia et al., 2013).

Weighted undirected matrices were analyzed using Brain Connectivity Toolbox (Rubinov and Sporns, 2010), igraph (version 1.5.1.9005; Csardi and Nepusz, 2006) in R, and tools in MATLAB (Mathworks, Natick, MA). Global graph metrics were calculated for edge density thresholds ranging from 6 to 32%. Global network measures for this density range were converted to area under the curve (AUC) values prior to multivariate statistical assessments. Node-specific (ROI) network measures were calculated at a 16% threshold, which in the present dataset corresponds to a Pearson r of approximately 0.2 or more, an average degree of 9.1, and network strength of 2.5–5.0 (Supplementary Figure 1). As a comparison, graph density based on bootstrapped calculations of confidence intervals for significant Pearson correlations corresponds to a significant r of ≥0.06, 40% density (δ ≥ 0.4), degree(k) ≥ 20, and strength ≥4 (Supplementary Figure 1).

A major objective of this study was to assess maturational and aging changes in functional network topology and determine links between network measures and spontaneous ‘affective’ behaviors. To this end, we assessed several graph measures of network integration and communication efficiency. These included node clustering coefficient (CC; number of interconnected node pairs connected to a node; Onnela et al., 2005), node strength (sum of edge weights/node) and degree (sum of edges/node) and global measures, such as transitivity (related to CC; number of triad groups normalized by all possible triad nodes in a network), mean distance or characteristic path length (CPL; the average edges or edge weights between node pairs), the small world coefficient (SWC; Humphries et al., 2006; Humphries and Gurney, 2008) and small world measure (SWM; Telesford et al., 2011). Global efficiency is the inverse of CPL. For local node efficiency, a length matrix was used to calculate inverse values in vicinity of a node, with added weights used to emphasize the highest efficiency node paths (Rubinov and Sporns, 2010; Wang et al., 2017). To corroborate results relative to random networks, all network measures were calculated on original and randomized versions of the functional connectivity matrices. Positive and negative edges were randomized by ~5 binary swaps with edge weights re-sorted at each step. The approach preserves weights, degree, and strength distributions (Rubinov and Sporns, 2010). We confirmed that the randomization approach preserved weights per nodes, node degree and strength distributions (Supplementary Figure 2). Node strength and degree distributions were similar between random and original networks, with distributions having long tail at low density threshold and near Gaussian at high density (Supplementary Figures 2A–F). Degree and strength were approximately preserved across nodes, although the specific organization of edge weights differed when viewed on 3D brains (Supplementary Figures 2G–K). The positive and negative matrix strength sequence accuracies for input and output matrices were r = 0.99 and r = 0.98. SWC was calculated as the ratio of gamma to lambda coefficients, where gamma is the ratio of real to random network CC, and lambda the ratio of real to random network CPL (Humphries and Gurney, 2008). SWM was calculated with an artificial non-ring lattice network as reference, as previously reported (Telesford et al., 2011).

A probabilistic approach for community detection was used to calculate a modularity statistic (Q), which indexes the rate of intra-group connections vs. connections due to chance (Blondel et al., 2008). The procedure starts with a random grouping of nodes and iteratively moving nodes into groups which maximize the value of Q. The final number of modules and node assignments to each group (e.g., community affiliation assignments) was taken as the median of 100 iterations of the modularity maximization procedure (Pompilus et al., 2020). The number of communities per mouse brain functional network and population size of each community were also assessed. We evaluated other community detection routines in igraph as a comparison to results obtained with the implementation of Louvain’s algorithm in Brain Connectivity Toolbox (results in Table 1). The role each node plays within their assigned module can vary according to their participation coefficient. The participation coefficient is calculated as: P=1−∑mimndimdi2, where d_im is the within-module node degree and d_i is the overall network degree for a node.

We analyzed the tendency of assortative vs. dissortative mixing of nodes (Newman, 2002). The assortativity index is a correlation coefficient comparing node strength values between pairs of edge-connected nodes. Positive r values indicate connectivity between pairs of nodes with similar strengths (e.g., high strength nodes pairs with high and low with low), while negative r values indicate cross-pairings between low and high strength nodes. We generated and analyzed weighted rich club coefficient curves (Φ) at an edge density of 16% (Colizza et al., 2006). Nodes belonging to a rich club subnetwork have an above-chance tendency to tightly connect with each other and share the bulk of connections within the network (Colizza et al., 2006; van den Heuvel and Sporns, 2011). The approach creates subgraphs containing nodes with strength values at a predetermined degree value, k. For each k-level subgraph, Φ_w is calculated as the ratio of the sum of edge weights that connect nodes of degree k or higher (W_>k) to the total possible number of summed edge-weights (van den Heuvel and Sporns, 2011). Rich club coefficient values were calculated for original and randomized networks and a normalized rich club index analyzed (i.e., Φnorm=Φorig/Φrand).

Functional connectivity networks were visualized in BrainNet viewer (Xia et al., 2013). Center coordinates for each node were based on CCFv3 parcellations, as indicated above. Several 3D mouse brain connectome maps were generated in which size of nodes (spheres) was scaled according to AUC values for several node-centric measures at various thresholds. Additional 3D maps were generated with node colors representing the module assignment (e.g., community affiliation vector) and node size weighted by modularity index.

Fresh brain tissue was extracted, and two isolated regions were weighed, and flash frozen for metabolomics (Supplementary Figure 3). A 1 mm thick slice of the prefrontal cortex was manually dissected from approximately Bregma anterior–posterior −2.96 mm to −1.94 mm (Paxinos and Franklin, 2001), and a 2 mm wide midline cut including from the median −/+ 1 mm medial lateral structures −3 to −3.5 mm deep. The dissected regions thus included areas of the infralimbic, prelimbic, medial orbital, anterior cingulate and midline primary motor cortices. The dorsal hippocampus was dissected at −1.46 mm rostrally to about −2.46 mm caudally (Bregma coordinates) and stripped of its overlying cortex and white matter. These regions have been implicated previously in age related differences in structural connectivity (Clifford et al., 2023).

Mice brain tissue (20–25 mg) was extracted with 1 mL of Acetonitrile:Isopropanol:Water (3:3:2, v:v:v), centrifuged for 25 min at 10,000× g. The supernatant was transferred into a new tube and dried in a speedvac. The dried extract was resuspended in 0.5 mL of Acetonitrile:Water (1:1, v:v) and centrifuged for 15 min and the supernatant was transferred into reaction vial and dried down under nitrogen. Fifty microliters of 2% MOX reagent were added to the vial and incubated at 30°C for 90 min followed by addition of 50 μL of MTBSTFA +TBDMS and incubation at 60°C for 1 h. Finally, reaction solution was transferred into GC vial for GC–MS analysis (Mahar et al., 2021).

GC–MS data was recorded with Single Quadrupole Mass Spectrometer and GC (Trace 1,310; Thermo Fisher Scientific, United States). The column was 30 m long dimethyl (95%)/diphenyl polysiloxane (5%) RTX-5MS with 0.25 mm ID, and guard column (10 m; Restek). The initial temperature of GC oven was 60°C for 60 s and the temperature increased to 325°C (@10°C/min). Metabolites were identified using NIST library with Xcalibur software.

Statistical analyses used MATLAB, GraphPad Prism (version 10.2) or the R statistics package. Unless otherwise noted, all analyses involved either a two-way full factorial analysis of variance (ANOVA) without or with repeated measures, depending on the specific test (critical α < 0.05). Post-hoc tests used Tukey’s honest significant difference (HSD) procedure. Where appropriate, false discovery rate (FDR, q ≤ 0.05) correction was used. The resultant statistical test parameters are indicated in figure legends or in the results section. Unless otherwise specified, male and female data were pooled to increase statistical power. The dependent variables in locomotor tests were sequential and repetitive beam breaks, total, margin, and center distances. Dependent variables for sucrose preference tests are fluid intake consumed in 24 h (ml) and sucrose preference index. For social interaction, the dependent variables are number of ROI entrances, average speed in ROI (in pixels per frame, or ppf), time spent in ROI, and time of initial approach toward ROI. For cued fear conditioning, the dependent variable is freezing counts summed over four 60 s trial bins per day.

Global network measures (i.e., SW, modularity, CPL, efficiency, transitivity, assortativity) were analyzed first as AUC values. Measures with significant differences were again analyzed across density thresholds to determine the threshold level where significant differences are observed. Upon determining the significance of modularity and assortativity, we analyzed the number of modules per group, number of nodes per module, and the normalized rich club index. All single threshold nodal analyses were carried out for a network density threshold of 16%. These included node degree, node strength, participation coefficient, clustering coefficient, and local efficiency. These were corrected for multiple comparisons across 60 nodes using the Benjamini-Hochberg FDR procedure (q ≤ 0.05).

Additional exploratory analyses were carried out on several variables to assess novel links between network measures and affective behavioral measures. For this, we used 9 variables, including AUC values for SW, modularity, CPL, transitivity, assortativity, freezing counts on days 1–3 and margin distance. This data set provided 9 variables with 30 observations each. Factor analysis was first used to explore associations between variables of most importance in the data set (i.e., evaluate variable loadings for the first two factors contributing most to the variance in the data). Projections of variables and data points onto the first 2 dimensions was carried out with principal components analysis (PCA) and cluster analysis using FactoMineR (Lê et al., 2008). A multivariate ANOVA (MANOVA) was used to examine overall differences and report specific univariate associations of potential future interest (α < 0.05).

Further analysis of fMRI scans was carried out using probabilistic ICA (Beckmann and Smith, 2004) as implemented in MELODIC Version 3.15, part of FSL. The following data pre-processing was applied to the input data: masking of non-brain voxels; voxel-wise de-meaning of the data; normalization of the voxel-wise variance. Pre-processed data were whitened and projected into a 20-dimensional subspace using PCA. The whitened observations were decomposed into sets of vectors which describe signal variation across the temporal domain (time-courses), the session/subject domain and across the spatial domain (maps) by optimizing for non-Gaussian spatial source distributions using a fixed-point iteration technique (Hyvarinen, 1999). Estimated component maps were divided by the standard deviation of the residual noise and thresholded (0.6) by fitting a mixture model to the histogram of intensity values (Beckmann and Smith, 2004). The resulting components were overlaid on the multi-subject template and classified according to peak z score anatomical location. A series of two-stage multiple linear regressions were used to back-project the group ICA components to subject-specific time courses (using spatial regressions) and spatial components (using temporal regressions). The subject-specific spatial component maps were then used in statistical comparisons for each component in the FSL randomize tool. The statistical design matrix was generated using the FSL Glm tool. We used a one-way ANOVA general linear model design. Randomization tests were carried out with 500 permutations (corrected p-level for significance is 0.05 and statistical thresholding of maps by threshold free cluster enhancement).

For metabolomics data, the intensity of metabolites was imported to the web-based MetaboAnalyst 5.0 platform¹ and normalized by the sum of intensities followed by generalized log10 transformation and pareto scaling before multivariate analysis. Further, principal component analysis (PCA), partial least square discriminant analysis (PLS-DA), and Hierarchical Clustering analysis were employed on the dataset (Pang et al., 2021). The Q² and R² values were calculated to assess supervised PLS-DA model. The Variable Importance in Projection (VIP) derived from PLS-DA was used to extract the importance of metabolites to PLS-DA model. Metabolites with a VIP score > 1 is weighed to be influential in PLS-DA.

Mice were tested in a novel open field environment 1 week after arrival. The ages at the time of testing were 11, 31, and 61wo. We used a one-way ANOVA with Tukey’s multiple comparison post-hoc test to analyze the data. The 10wo group showed a significantly greater total distance (F_2,27 = 4.5, p = 0.02; Tukey’s post hoc p = 0.03 and 0.04 comparing 10wo vs. 30–60wo, respectively) and margin distance than 30 and 60wo mice (F_2,27 = 6.3, p = 0.006; Tukey’s post hoc p = 0.01 comparing 10wo vs. both 30–60wo mice; Figure 3). The greater margin distance activity in 10wo mice was driven by female mice (sex main effect F_1,24 = 5.9, p = 0.02, with margin distance higher in 10wo female mice compared to 30–60wo females and 10wo males; Tukey’s post hoc p < 0.05; Supplementary Figure 4). No significant age or sex differences were observed in other locomotor activity measures, such as sequential and repetitive beam breaks, total and center distance.

Sociability tests were carried out when mice were 12, 32, and 62wo. We used a two-way ANOVA (age x test session) with repeated measures. No effect of sex was observed. Thus, data were pooled in our analyses. No age differences in the number of ROI entrances was observed but a significant difference across test sessions was observed (main effect of session F_2,46 = 49.4, p < 0.00001; Figure 4A). Thus, all age groups entered the ROI where the novel test mouse was located a similar number of times and displayed reductions in the number of entrances during a second encounter with the same test mouse by the same degree. 10wo mice showed a greater mean speed within the ROI than 60wo mice, during the re-encounter with the familiar mouse and exploration of the 2nd novel mouse (F_2,23 = 5.9, p = 0.008; Tukey’s post hoc p = 0.02 and 0.03, respectively; Figure 4B). We also observed significant effects of test session and age, but not age x session interactions, on time spent in the test mouse ROI (session F_2,46 = 5.5, p = 0.007; age F_2,23 = 4.5, p = 0.02; Figure 4C). 30wo mice spent more time in the ROI during the first novel mouse encounter than 10wo mice (Tukey’s post hoc p = 0.007). Only 30wo mice showed reductions in their time spent in the ROI on the second encounter with the familiar mouse and with the new novel mouse (Tukey’s post hoc p = 0.02 and 0.005, respectively). No effects of age or session were observed on time in initial approach (Figure 4C).

Sucrose preference testing was conducted a week after sociability tests when mice were 13, 33, and 63wo. No effect of sex was observed. Thus, data were pooled in our analyses. All mice consumed more sucrose solution than water. No age-related differences in the amount of sucrose and water consumed nor differences in the sucrose preference index were found (Supplementary Figure 5). We also analyzed sucrose consumption during the 3 days of acclimation and observed no differences between the groups (Supplementary Figure 5).

Fear conditioning was conducted a week after completing fMRI scans when mice were 16, 36, and 66wo. Percent time freezing was analyzed using a two-way repeated measures ANOVA. No effect of sex was observed. Thus, data were pooled in our analyses. On day 1 of fear conditioning, all age groups progressively developed unconditioned-conditioned stimulus (US+/CS+) associations (Figure 5). We observed a rise in percent time freezing over the course of each tone-shock exposure in all age groups (effect of trial F_3,81 = 66, p < 0.0001 and age F_2,27 = 7.5, p = 0.0025; Figure 5A). 60wo mice showed higher percent freezing across cue-shock trials than the other age groups (Tukey’s post-hoc test indicated p < 0.05 comparing 60wo to 10wo and 30wo mice on trials 1–3 and p = 0.02 compared to 10wo mice on trial 4), although the freezing responses were more variable in the 60wo mice than in the other two age groups. Data were next analyzed as average percent time freezing per day (Figure 5B). When compared across days, all age groups showed robust recall of fear conditioning, observed as high levels of freezing in response to the tone alone, at 24 and 48 h (main effect of age F_2,27 = 4.6, p = 0.02; main effect test day F_2,54 = 66.8, p < 0.0001). At 24 h, all mice maintained a high percent freezing relative to day 1 (Dunnett’s test p < 0.05 in all age groups comparing day 2 to 1). At 48 h, the cages were modified such that the tone was delivered in a context with different sensory features (walls and floors were smooth white plastic in a triangular configuration, as reported in Anagnostaras et al., 2010). Still, 10wo and 60wo mice maintained high levels of freezing in response to the conditioned cue, relative to day 1 (Dunnett’s test p < 0.05). The 30wo mice showed a similar but non-significant trend toward increased fear recall at 48 h (Figure 5B). The freezing responses by 60wo mice were higher than those of 10wo mice on the 3^rd day of testing (Dunnett’s test p < 0.05).

We analyzed a motion index as a surrogate measure of general locomotor activity using an infrared motion sensitive camera inside the test cages. Consistent with results in Figure 3, 10wo mice showed higher motion index levels than 60wo mice in the novel test cage on day 1 (Figure 5C; day by age interaction F_4,54 = 2.7, p = 0.04, Tukey’s HSD p = 0.008 day 1 comparing 10wo vs. the 60wo group; main effect test day F_2,54 = 84.9, p < 0.0001). The motion index was significantly reduced in all groups by days 2 and 3 of testing (Tukey’s HSD p < 0.05 comparing days 1 vs. 2 and 3 for all age groups). A closer assessment of the motion index over the tone-shock presentation trials on day 1 and tone responses on day 2 provided further evidence of highly variable responses to conditioned and unconditioned stimuli in 60wo mice compared to the other two age groups (Supplementary Figures 6–8). The 60wo mice did have a reduced initial motion response to the tone, but this was only significant during the second tone-shock presentation (two-way ANOVA with repeated measures showed effect of time, F_2.8,76.1 = 46, p < 0.0001, and age, F_2,27 = 7.5, p = 0.0025, but no interactions; Tukey’s multiple comparison test indicated p < 0.05 comparing 60wo to 30wo, on tone-shock trial 2). However, all mice responded to tone, shock, and showed reduced motion during baseline over the course of all 4 trials on day 1 (Supplementary Figures 6D, 8). Further declines in motion were not observed on day 2 (Supplementary Figure 7D), suggesting acclimation of motion overlapping with conditioned fear learning.

Functional MRI scanning was carried out over a two-week period following sucrose preference testing. Mice were 14–15, 34–35, and 64–65wo at the time of scanning. Data were analyzed using two-way ANOVA with repeated measures and Tukey’s multiple comparison test. No effect of sex was observed. Thus, data were pooled in our analyses. Analysis of network measures obtained from functional connectivity matrices indicated significant main effects of age, but not sex, for SWC (age by threshold interaction F_26,312 = 2.4, p = 0.0002), modularity (age by threshold interaction F_26,312 = 1.79, p = 0.01), and assortativity (age by threshold interaction F_26,312 = 2.6, p = 0.00006; Figure 6). Tukey’s post-hoc comparisons further revealed that SWC was highest in 30wo mice relative to 10wo mice and both modularity and assortativity were highest in 10wo mice relative to 30 and 60wo mice (Figure 6). Significant pairwise differences were observed across several edge density thresholds (see Figures 6C–E; details in figure legend). We analyzed the number of modules (at 16% edge density) per age group and observed a significantly greater number of modules in 30 vs. 10wo mice (Figures 6A,F; Kruskall-Wallis test p = 0.02). Network randomization removed these differences (Figure 6B). The number of nodes per module did not differ between groups (Figure 6G). To confirm the significant age related changes in functional network modularity, we further investigated 5 additional community detection routines available in igraph, including greedy optimization (Clauset et al., 2004), edge betweenness (Newman and Girvan, 2004), random walk (Pons and Latapy, 2005), lead Eigenvector (Newman, 2006), and a different implementation of Louvain’s algorithm (Blondel et al., 2008; Table 1). These were all estimated at 16% edge density. Across most of these algorithms, the number of modules detected was lowest in 10wo mice compared to the other two age groups. In addition, the modularity index was highest in 10wo mice compared to the two older age groups (Table 1 provides statistical results).

At 16% edge density, we also observed significant main effects of age on transitivity and assortativity (Table 1). Transitivity was lowest in 60wo vs. 30 and 10wo mice. Assortativity was significantly greater in 10wo vs. 30wo mice. This was observed using Brain Connectivity Toolbox’s routine for calculating network assortativity. Although a similar trend was observed using igraph, it did not reach statistical significance (Table 1).

Group ICA identified 9 of 20 functional network components similar to previously identified in mice (Jonckers et al., 2011; Sforazzini et al., 2014; Grandjean et al., 2014a; Bukhari et al., 2017). The ICA components included anterior and posterior cingulate, somatosensory, visual, auditory/temporal, piriform/amygdala, hypothalamic and dorsal striatal networks (Supplementary Figure 9). Figure 7 shows the results for statistical comparisons between age groups. The 10wo mice had greater functional connectivity in the auditory/temporal network than 60wo mice. The areas showing greater connectivity within the right hemisphere of this network in 10wo mice included anterior cingulate, auditory/temporal subregions, ventral hippocampus, dorsal striatum, midbrain reticular nucleus, bed nucleus of stria terminalis, posterior thalamic areas, and superior colliculus. The left hemisphere had less differences and included primary and secondary somatosensory areas, striatum, visual cortex, hippocampus, and thalamic subregions (Figure 7). The left hemisphere striatal network had greater functional connectivity in 10wo compared to 60wo mice. The greater functional connectivity in the left striatal network of 10wo mice included piriform, hippocampus and superior colliculus. We also observed greater functional connectivity in a left hemisphere striatal network in 30wo compared to 60wo mice and in a right hemisphere somatosensory network in 10wo compared to 30wo mice.

MANOVA revealed a significant effect of age across several connectomic and behavior variables (F_9,20 = 4.6, p = 0.002). Principal components analysis (PCA) revealed a relatively strong relationship between several of these variables. This is observed by contributions to the first 2 dimensions, which jointly explain 56.08% of variance in the data. To establish significance, PCA results for the first 2 dimensions (56.08%) were compared against a reference value (43.65%) representing 95% quartiles of distribution of percent explained variances (inertia) of simulated data tables of equivalent size and normal distribution (Figure 8). The positive direction along dimension 1 opposes 90% of mice in the youngest group (1–8 and 10 are in the youngest 10mo group) from many of the mice in the other two age groups (Figure 8A; Supplementary Figure 10A). This positive direction is characterized by high values freezing counts on days 2–3 of fear conditioning and high margin distance activity in a novel open field environment and high values for small world coefficient, transitivity, and modularity in 70% of the 10mo mice (Supplementary Figure 10B). Cluster dendrogram in Figure 8B classified 3 clusters. A majority of young 10mo mice were classified in cluster 3 (green). This group is again characterized by high values for modularity, transitivity, small world coefficient, margin distance activity and days 2–3 of fear conditioning.

Results for PCA and partial least squares discriminant analysis (PLS-DA) for GC–MS data obtained from the mouse hippocampus are shown in Supplementary Figure 11 and Figures 9A,B. Variable importance in projection (VIP) scores plot of the PLS-DA model for mouse hippocampus indicated differences in the top 25 metabolites between 10, 30, and 60wo mice. Among the top 25 listed in Figure 9A, methionine and cysteine had the highest VIP scores. The top 7 had a bidirectional transition in relative levels when youngest mice having lowest or highest levels and late middle-aged mice having diametrically distinct levels. This relationship included methionine, cysteine, malate, fumarate, succinate, N-acetylglycine, alanine, and taurine. The differential prominence of these metabolites across age groups is supported by clustering in heatmap shown in Figure 9B. Results for PCA and PLS-DA for GC–MS data obtained from the mouse midline prefrontal cortex are shown in Supplementary Figure 12 and Figures 9C,D. VIP scores plot of the PLS-DA model for mouse prefrontal indicated differences in the top 25 metabolites between 10, 30 and 60wo mice. Among the top 25 listed in Figure 9C, proline, taurine, and methionine had the highest VIP scores. The differential prominence of these metabolites across age groups is supported by clustering in the heatmap in Figure 9D.