Animals
All procedures involving mice were performed according to protocols approved by the Stanford University IACUC and APLAC (protocol no. 8661) and VA Palo Alto Committee on Animal Research ACORP (LUO1736). For the ageing cohort and exercise cohort, male C57BL/6JN mice were obtained from the National Institute on Aging (NIA) Aged Rodent colony. For the whole-body partial reprogramming cohort, male whole-body inducible OSKM (iOSKM) mice (ROSA26(rtTA-M2); Col1a1(tetO-OSKM)) (on a mixed background of the following strains: C57BL/6, B6D2F1, 129S4 and B6129SF1/J) were generated from the Jaenisch laboratory80 and obtained from the Jackson Laboratory (JAX 011004). Mice were housed in groups of 3–5 mice of the same age at the ChEM-H/Neuro vivarium (ageing and partial reprogramming) or at the Veterinary Medical Unit at the Veterans Affairs Palo Alto Health Care System (exercise) under 12-h light/dark cycles, approximately 21 °C ambient temperature, and approximately 50% humidity for at least 3 weeks before any experiments or sample collection occurred.
Ageing coronal and sagittal cohorts
We used two independent ageing cohorts of male C57BL/6JN mice (referred to as ageing cohort 1 and ageing cohort 2). Ageing cohort 1 included mice of the following ages for the coronal section dataset: 3.8, 5.4, 9.8, 12.9, 15.5, 21.4, 23.5, 26.7, 30.9 and 33.2 months. Ageing cohort 1 also included mice of the following ages for the sagittal section dataset: 3.8 and 26.7 months. Ageing cohort 2 included mice of the following ages for the coronal section dataset: 3.4, 4.3, 6.6, 15.8, 18.8, 19.8, 24.6, 28.5, 32.6 and 34.5 months. Ageing cohort 2 included mice of the following ages for the sagittal section dataset: 6.6, 8.6, 19.8 and 23.5 months. Information on all mice in the ageing cohorts is provided in Supplementary Table 4.
Exercise experiment
The exercise experiment included 3 groups of male C57BL/6JN mice: 4 young (3 months) sedentary mice, 4 old (19 months) sedentary mice, and 4 old (19 months) exercise mice. Information on all mice in the exercise experiment is provided in Supplementary Table 4. Sample sizes were selected to allow testing of statistically significant differences of mouse-level attributes (for example, cell-type proportion) using the non-parametric two-sided Mann–Whitney U-test across groups. Littermates were randomized for exercise and sedentary conditions. The mice in the old exercise group were provided with voluntary wheel running through individual housing for 5 weeks in polycarbonate cages with 12.7 cm diameter running wheels (Lafayette Instrument, 80820) and monitored weekly for adequate running. The five-week duration of voluntary exercise was selected based on published regimens for efficacy32,81. Sedentary mice were individually housed in identical ages without running wheels.
Whole-body partial reprogramming experiment
The whole-body partial reprogramming experiment included 3 groups of male iOSKM mice (on a mixed background of the following strains: C57BL/6, B6D2F1, 129S4 and B6129SF1/J): 4 young (4.8–4.9 months) control mice, 4 old (25.6–29.2 months) control mice, and 4 old (26.5–29.2 months) OSKM mice. Information on all mice in the partial reprogramming experiment is provided in Supplementary Table 4. Sample sizes were selected to allow testing of statistically significant differences of mouse-level attributes (for example, cell-type proportion) using the non-parametric two-sided Mann–Whitney U-test across groups. Old mice were matched by age and body weight and then randomized to control and OSKM conditions. All mice (control and OSKM) were individually housed during the experiment. The mice in the old OSKM group underwent three periods of cyclic induction of OSKM by doxycycline treatment, which consisted of doxycycline administration in the drinking water for two days (ON), followed by five days without doxycycline administration (OFF) and repeated for 3 weeks (ON–OFF–ON–OFF–ON) with mice euthanized at the end of the last doxycycline administration treatment. Doxycycline (Fisher ICN19895505) was dissolved in drinking water (1 mg ml−1), placed in amber water bottles to protect the solution from light, and provided ad libitum to the mice in the old OSKM group.
Sample collection
Mice were euthanized to collect fresh frozen whole-brain samples for MERFISH experiments. Mice were euthanized with 5 min of exposure in a CO2 chamber, and brains were removed and placed in a cryomold on ice and filled with pre-chilled Optimal Cutting Temperature (OCT) Compound (Fisher Healthcare Tissue Plus, 4585) and then placed on dry ice. After OCT solidified, the samples were moved to long-term storage at −80 °C. RNaseZap (Invitrogen, AM9780) was used to disinfect all dissection tools before and after each mouse. The sample collection for ageing cohort 1 occurred on 4 November 2022 from 14:00–15:30 pm Pacific Standard Time (PST) and the sample collection for ageing cohort 2 occurred on 28 June 2023 from 14:30–16:15 PST. Sample collection alternated between younger and older mice. For the exercise experiment, mice were perfused with 15 ml of PBS with heparin sodium salt (50 U ml−1) (Sigma-Aldrich, H3149-50KU) before sample collection. Sample collection for the exercise experiment occurred on 6 June 2022, and sample collection for the partial reprogramming experiment occurred on 28 June 2023. Sample collection alternated between mice from different experimental conditions.
MERFISH 300-gene panel selection
We selected 300 genes to profile for the MERFISH experiments. Our selected 300-gene panel consists of 129 cell-type and subtype markers (81 cell-type markers, ranging from 1 to 8 markers per cell type; 48 subtype or function-related markers) as well as 181 genes of interest that have been implicated in ageing-related pathways, in important processes not previously linked to ageing, as well as genes identified by analysis of single-cell RNA-seq atlases of murine brain ageing, with 10 genes shared across cell-type and subtype markers and genes of interest (Supplementary Table 1).
The cell-type and subtype marker genes included all markers in a suggested panel by Vizgen; markers determined from literature review focused on cells of the brain vasculature5,82,83, NSCs and neuroblasts3,5,84,85, and immune cells3,5,86,87. For these 129 cell-type or subtype markers, we also included several NSC and neuroblast markers obtained from integrated differential gene expression analysis across multiple single-cell RNA-seq datasets of the adult mouse subventricular zone3,4. We included markers for the following cell types: excitatory neurons, inhibitory neurons, medium spiny neurons, astrocytes, microglia, oligodendrocytes, OPCs, endothelial cells, pericytes, VSMCs, VLMCs, ependymal cells, neuroblasts, NSCs, macrophages, neutrophils, T cells, B cells, natural killer (NK) cells, mast cells and dendritic cells.
For the ageing-related genes, we included a set of 33 genes related to murine ageing selected from the GenAge model organism database28 (accessed 7 September 2022) (Supplementary Table 1) to enrich for genes known to have a causal role in ageing, which is helpful for building more causal ageing clocks88,89. We further included genes related to T cell activity3, subventricular zone NSC heterogeneity90, endothelial heterogeneity91, meningeal lymphatics functions83,92,93, cellular senescence, immune response, stem cells, neurogenesis94 and vasculogenesis95. We also included sets of expert-curated genes pertaining to interesting cellular and organismal functions, some of which have emerging roles in the regulation of ageing but are not well studied yet, including T cell signalling, reprogramming, cell adhesion and migration, lipid metabolism and neuropeptide signalling. Finally, for a more unbiased set of genes, we included several differentially expressed genes (DEGs) between young and old mice for cell types across three single-cell transcriptomics datasets of the subventricular zone3,4 and two multi-region brain single-cell transcriptomics datasets5,7.
We ensured that all selected genes in the MERFISH panel were expressed in existing single-cell or single-nuclei RNA-seq brain atlases4,5,7 and met technical constraints for minimizing optical clouding. This included limiting the total estimated gene expression using the Vizgen Gene Panel Design Portal to under 9,000 fragments per kilobase of transcript per million mapped reads (FPKM) (total estimate at 7,691 FPKM) and limiting the maximum estimated expression per gene using the Vizgen Gene Panel Design Portal to under 700 FPKM (maximum per-gene estimate at 452 FPKM). The complete gene panel, classification of markers, and rationale for inclusion are included in Supplementary Table 1.
MERFISH imaging experiment
The MERFISH experiment was conducted through the Vizgen MERSCOPE technology laboratory service. Compared with other spatially resolved single-cell transcriptomics, Vizgen MERSCOPE technology has been shown to provide high specificity and sensitivity even with larger gene panel sizes96. Fresh frozen mouse brain samples were cut into 10 μm-thick sections on a cryostat at −20 °C and placed onto a MERSCOPE slide (Vizgen 20400001). We obtain coronal sections that contain the CTX, STR, CC/ACO and VEN including the subventricular zone neurogenic niche. We obtain sagittal sections that contain the aforementioned brain regions along with additional brain regions (olfactory bulb, rostral migratory stream, brain stem, and cerebellum). Two coronal sections were placed on each slide and paired to balance ages while sagittal sections were placed on their own slides. The tissue sections were fixed with 4% paraformaldehyde in 1× PBS for 15 min, washed 3 times with 5 ml 1× PBS and incubated with 70% ethanol at 4 °C overnight for tissue permeabilization. Samples were then stained for cell boundary using Cell Boundary Kit (Vizgen, 10400009), and later hybridized with a custom designed MERSCOPE Gene Panel Mix consisting of 300 genes (Vizgen 20300008) in a 37 °C incubator for 36–48 h. Following incubation, the tissues were washed with 5 ml formamide wash buffer at 47 °C for 30 min, twice and embedded into a hydrogel using the Gel Embedding Premix (Vizgen 20300004), ammonium persulfate (Sigma, 09913-100 G) and TEMED (N,N,N′,N′-tetramethylethylenediamine) (Sigma, T7024-25ML) from the MERSCOPE Sample Prep Kit (10400012). After the gel mix solution solidified, the samples were cleared with clearing solution consisting of 50 μl Proteinase K (NEB, P8107S) and 5 ml of Clearing Premix (Vizgen 20300003) at 37 °C overnight. After removing clearing solution, the sample was stained with DAPI and Poly T Reagent (Vizgen 20300021) for 15 min at room temperature, washed for 10 min with 5 ml of Formamide Wash Buffer, and then imaged on the MERSCOPE system (Vizgen 10000001). A fully detailed, step-by-step instruction on the MERFISH sample prep the full protocol is available at https://vizgen.com/resources/fresh-and-fixed-frozen-tissue-sample-preparation/. Full instrumentation protocol is available at https://vizgen.com/resources/merscope-instrument/. The MERFISH data were collected over two separate batches (see ‘slide_id’ in Supplementary Table 4). The first batch, A, included ageing cohort 1 (both coronal and sagittal) and half of the exercise experiment samples with even representation across conditions. The second batch, B, included ageing cohort 2 (both coronal and sagittal), the remaining half of the exercise experiment samples, and the partial reprogramming experiment samples.
Cell segmentation and MERFISH data preprocessing
Segmentation of cells was performed using Cellpose (1.0.2) through Vizgen’s laboratory service. Cell segmentation was implemented on images using nuclear staining (DAPI) and cytosolic staining (Poly T). Transcripts were allocated using these cell segmentations by summing across seven z-stacks, accounting for both nuclear and cytosolic (soma) transcripts. Quality control statistics were computed using the Vizgen post-processing tool (vpt) (1.2.2) for cell segmentation (Supplementary Table 2).
For preprocessing, we performed initial cell filtering separately for each MERFISH experiment (consisting of either two coronal sections on the same slide or one sagittal section). For each experiment, we removed putative doublets using Scrublet97 and a doublet score cutoff of 0.18. We then filtered out all cells with segmentation volume less than or equal to 100 μm3 or greater than or equal to 3 times the median cell volume. We also filtered out all cells with fewer than or equal to 20 counts and/or fewer than or equal to 5 genes with non-zero expression. To correct for potential different segmentation sizes, we divided the raw transcript counts obtained for each cell in the MERFISH dataset by the volume of the corresponding segmentation. After combining all experiments into an integrated dataset, we then filtered out all cells in the top 2% highest and top 2% lowest total expression. Statistics associated for the aforementioned cell filtering steps and for additional steps after clustering can be found in Supplementary Table 3. To obtain log-normalized gene expression values, we normalized the total gene expression for each cell to 250 and log-transformed the expression with an added pseudocount. This procedure was performed separately for the ageing cohorts (coronal), ageing cohorts (sagittal), exercise experiment and partial reprogramming experiment.
Cell-type clustering and identification
For clustering, we converted each log-normalized gene expression value to a z-score using scanpy.pp.scale with max_value = 10 in the scanpy package98. We performed Leiden clustering using scanpy.tl.leiden with resolution = 0.5 for the initial clustering and default settings otherwise. We obtained batch-balanced nearest neighbours graph using BBKNN (scanpy.external.pp.bbknn) and then used this neighbours graph to generate a UMAP visualization of all cells (scanpy.tl.umap). For the partial reprogramming experiment, which only involved one batch of MERFISH data, we used scanpy.pp.neighbors with n_pcs=20 and n_neighbors=15 instead of BBKNN. To annotate cell types, we manually labelled each cluster based on cell-type expression patterns as observed across two orthogonal data visualization modalities (the UMAP visualization and a heat map of cell-type markers) to reduce errors in cell-type annotation resulting from dimensionality reduction distortions99,100. For clusters that expressed markers from multiple cell types, we performed successive Leiden clustering on those clusters until unique cell types could be annotated (see Supplementary Table 5). This procedure was performed separately for the ageing cohorts (coronal), ageing cohorts (sagittal), exercise experiment, and partial reprogramming experiment. A detailed description of the cell-type markers and Leiden clustering resolutions for each dataset and cell type can be found in Supplementary Table 5. Although a small set of markers for some rare immune cell types (NK cells, mast cells, dendritic cells) were included in our MERFISH panel (Supplementary Table 1), we were unable to identify these cell types in any of the datasets (Supplementary Table 5), probably owing to their low abundance. NK cells, mast cells, and dendritic cells were also not identified in several existing spatial transcriptomics studies of adult mouse brain17,21,22,23,24. NK cells and dendritic cells were identified in a large-scale and high-resolution spatial transcriptomic profiling of the whole adult mouse brain25, although at 3–4 times lower abundance than the rarest immune cell types identified in our datasets (T cells and B cells), which may explain why we were not able to identify these cells in our dataset. In addition, in the partial reprogramming dataset, we were unable to identify other rare immune cell types (T cells, B cells and neutrophils), consistent with our previously published dissociated single-cell RNA-seq datasets35 and probably owing to the lower abundance of these cell types in this partial reprogramming mouse model.
Spatial region and subregion clustering and annotation
To identify anatomical regions across the MERFISH datasets and assign region and subregion labels to each cell, we adapted a semi-supervised approach for clustering and annotating region labels from cell-type composition of local neighbourhoods around each cell17. For a given cell, we computed the cell-type abundances for each cell within 100 μm distance from the given cell. Then, we performed principal component analysis on the matrix consisting of the cell-type abundance profiles for each cell, applied k-means clustering (k = 25), manually visualized and merged clusters to obtain seven subregion annotations (CC/ACO, CTX_L1/MEN, CTX_L2/3, CTX_L4/5/6, STR_CP/ACB, STR_LS/NDB and VEN), and finally merged subregion annotations to obtain four region annotations (CC/ACO, CTX, STR and VEN). While we observed some variability in the subregion annotations between samples, there was general consistency in the four region annotations. We verified the expression of cortical layer markers in the three subregions of the CTX (Extended Data Fig. 3c,d), but we were unable to annotate each of the six known cortical layer individually through the clustering procedure, perhaps owing to the low number of cortical layer markers. This region and subregion clustering and annotation procedure was performed separately for the ageing cohorts (coronal), ageing cohorts (sagittal), exercise experiment, and partial reprogramming experiment.
Cell-type composition analysis
We computed cell-type proportions for a given sample by dividing the number of cells of each cell type by the total number of cells in the sample. For regional cell-type proportions, we divided the number of cells of each cell type in that region by the total number of cells in that region. Pearson correlation, 95% confidence interval for the correlation, and P value for association between cell-type proportion and sample age was computed using scipy.stats.pearsonr. We annotated strong changes in cell-type proportion with age as cell types with 95% confidence interval for the correlation that does not overlap with the interval [−0.25, 0.25]. Strong increases in cell-type proportion with age (coloured red) had 95% confidence interval for the correlation with lower bound greater than 0.25, and strong decreases in cell-type proportion with age (coloured blue) had 95% confidence interval for the correlation with upper bound less than −0.25. Linear regression of cell-type proportion on sample age with 95% confidence interval was computed using seaborn.regplot. To compute statistical significance of differences in cell-type proportions across categorical conditions, we used the two-sided Mann–Whitney U-test.
Increasing and decreasing gene expression with age analysis
For a given cell type, to identify genes that changed in expression with age, we computed the Spearman correlation between age and pseudobulk gene expression across samples in the coronal section dataset. The pseudobulk gene expression was computed as the mean log-normalized gene expression across all cells of the same cell type within a sample. For each gene, we obtained the Spearman correlation, the associated P value, and the lower and upper bounds of a 95% confidence interval for the correlation. We classified genes as ‘increasing’ if they had Spearman correlation greater than 0.3 and with the lower bound of the 95% confidence interval greater than 0.0. We classified genes as ‘decreasing’, if they had Spearman correlation less than −0.3 and with the upper bound of the 95% confidence interval less than 0.0. To reduce false positives resulting from transcript spillover due to segmentation, we constrained our analysis to the 220 genes with less than 5% estimated spillover based on an internal Vizgen metric using gene expression variations influenced by local cellular composition (see Supplementary Table 7 for list of genes).
GO enrichment analysis
We performed GO enrichment analysis to determine biological processes that were enriched in different sets of genes. For genes that increase or decrease in expression with age that were identified in ‘Increasing and decreasing gene expression with age analysis’, we performed GO biological process enrichment analysis separately for each set of genes and for each cell type. For genes in different spatiotemporal gene expression trajectory clusters for oligodendrocytes in the CC/ACO region (see ‘Spatiotemporal gene expression trajectory analysis’), we performed GO biological process enrichment analysis for all genes present in each of the nine trajectory clusters separately after filtering out genes with greater than 5% estimated spillover based on an internal Vizgen metric using gene expression variations influenced by local cellular composition. For GO enrichment analysis on genes used by the spatial ageing clocks, we selected positive coefficient clock genes (up to 50 genes with the largest positive coefficients) and negative coefficient clocks genes (up to 50 genes with the largest negative coefficients) and performed GO biological process enrichment analysis separately for each set of genes. For GO enrichment analysis on DEGs in endothelial cells in response to exercise in old mice, we selected genes that significantly increased with exercise (increased in old exercise compared with old sedentary with P < 0.05 from two-sided Mann–Whitney U-test) and genes that significantly decreased with exercise (decreased in old exercise compared with old sedentary with P < 0.05 from two-sided Mann–Whitney U-test) and performed GO biological process enrichment analysis separately for each set of genes.
We performed GO biological process enrichment analysis by selecting genes for each analysis (as described above), and using all other genes measured by MERFISH as background. GO enrichment was performed with topGO101 (R package version 2.54.0) using Fisher’s exact test for all Biological Process terms.
Regional gene expression changes with age
To compare different anatomic regions and subregions by the magnitude of gene expression changes with age, we selected the five youngest and five oldest mice in the data. For each cell type, we determined the minimum number of cells present across each of the mice and regions, excluding the VEN due to their low cell number. We then downsampled cells for each mouse and region, without replacement, to that minimum number, such that after sampling, all combinations of mouse and region had the same number of cells. We then excluded cell types in which the number of cells per mouse and region was less than 20. We calculated the transcriptional profile for each mouse and region by averaging across the volume-normalized expression of all cells from that mouse in that region, normalizing this profile to sum to 250, and performing a log transformation with an added pseudocount. To determine the change between old and young mice in a specific region, we subtracted the mean profile of the 5 young mice from the mean profile of the 5 old mice and averaged the absolute value of this difference across all genes after filtering out genes with greater than 5% estimated spillover based on an internal Vizgen metric using gene expression variations influenced by local cellular composition. We repeated this process 20 times, each time sampling different cells.
Spatiotemporal gene expression trajectory analysis
For each cell, we divided the raw transcript counts by the segmentation volume of the cell and then normalized by total count using the default settings of scanpy.pp.normalize_total. For each combination of cell type, anatomic subregion, and gene; we computed a vector of length 20 (trajectory) containing the mean expression of that gene in each age conditioned on a specified subregion and cell type. If the given cell type and subregion combination was not present in at least 70% of profiled ages, we classified that trajectory as ‘missing’. We performed last-value-carried-forward imputation to fill in missing values within trajectories that passed this threshold. Each trajectory was standardized by centring and scaling to unit variance with sklearn.preprocessing.StandardScaler102. Then, we performed k-means clustering using sklearn.cluster.KMeans with n_clusters=9, random_state=444, n_init=’auto’, and the matrix with each row corresponding to a scaled gene expression trajectory as input. The parameters of the clustering were selected to maximize the number of clusters while also maintaining qualitatively distinct trends in each cluster (Extended Data Fig. 4a,b). For each cluster of gene expression trajectories, we visualized the smoothed median and interquartile range of the scaled expression values. Smoothing was done with interpolating B-splines using scipy.interpolate.BSpline with s = 20. We manually annotated each cluster based on the qualitative expression patterns. These trajectory clusters include genes with expression that increases late in life (‘increasing late’), increases gradually throughout life (‘increasing gradual’), has the lowest expression in midlife (‘midlife trough’), peaks in early life (‘early peak’), peaks in midlife (‘midlife peak’), peaks in late life (‘late peak’), decreases after midlife (‘midlife decrease’), decreases early in life (‘decreasing early’), and decreases gradually through life (‘decreasing gradual’). To visualize representative gene expression trajectories for each cell type, we performed the same smoothing procedure separately for each cell-type-specific subset of trajectories from the clusters. We developed this trajectory clustering approach instead of using parametric methods (for example, polynomial fitting) to provide a more unbiased characterization of different trajectories.
SpatialSmooth soft pseudobulking procedure
For each cell type, we build spatial graphs connecting each cell with its 20 nearest neighbours by Euclidean distance and of the same cell type103. Spatial ageing clock performance was generally robust to the choice of the number of nearest neighbours (k) in building this graph (Extended Data Fig. 5a). We computed a L1-normalized adjacency matrix representing the spatial graph. Then, the SpatialSmooth algorithm propagates gene expression features across the cell-type-specific spatial graph by iterating the update equation until convergence (that is, \({X}_{t+1}\approx {X}_{t}\)):
$${X}_{t+1}\leftarrow (1-\alpha )X+\alpha S{X}_{t}$$
Where \(X\) is the initial gene expression matrix (cells as rows, genes as columns), \(S\) is the normalized adjacency matrix, and \(\alpha \) is the smoothing parameter. We set \(\alpha \) equal to 0.8. We set the convergence to be 30 iterations at maximum with a tolerance of 0.01, for which convergence is reached if \({| | {X}_{t+1}-{X}_{t}| | }_{\infty }\) is less than the tolerance or 30 iterations has elapsed. For rare cell types (for example, T cells, B cells, neutrophils), SpatialSmooth generally converged before 30 iterations. For common cell types (for example, oligodendrocytes, microglia, astrocytes), SpatialSmooth was generally performed for 30 iterations. After convergence at step \(t=T\), we use the smoothed spatial gene expression matrix \({X}_{T}\) as input for training the ageing clocks.
Training and cross-validated evaluation of spatial ageing clocks
For training spatial ageing clocks, we performed the SpatialSmooth procedure on the log-normalized gene expression to obtain smoothed spatial gene expression matrices for each cell type independently. Then, for each cell type, we fitted a pipeline consisting of standardization of gene features followed by lasso regression model to predict sample age from a cell’s gene expression profile and used sklearn.linear_model.LassoCV to select optimal hyperparameters with cv=5, n_alphas=20, max_iter=10000. We refer to this entire pipeline from SpatialSmooth to age prediction as the ‘spatial ageing clock’. To avoid explicit conditioning of age prediction on spatial information, the spatial ageing clocks only leverage the spatial context to process the input data via SpatialSmooth.
For cross-validated evaluation of spatial ageing clocks, we held out a single sample/age as the test set and kept the remaining samples/ages as the train set. SpatialSmooth was performed separately for the train and test sets. We fitted the lasso regression pipeline to predict age on the train set and used the model to obtain predicted ages on the test set. We repeated this procedure across all samples/ages to obtain predicted ages for all cells in the study. Performance was evaluated by Pearson’s correlation (R) and mean absolute error between the age and predicted ages of individual cells obtained from cross-validation. We also computed the Pearson’s correlation (r) between the age and median predicted ages of individual mice obtained from cross-validation. For the subregion-specific ageing clocks, we trained and evaluated the models using the same settings except with five nearest neighbours for SpatialSmooth and restricted to only cells in each subregion for training. To compare spatial ageing clocks to cell-type-specific ageing clocks trained on dissociated single-cell RNA-seq data from adult mouse subventricular zone, we used the median predicted ages associated with these clocks and data4.
Visualization of spatial ageing clock predictions
We used two approaches for visualizing the predicted ages obtained from the application of spatial ageing clocks, either through cross-validation on the coronal section dataset or directly through validation on an external dataset.
For datasets with relatively uniform distribution of many ages across lifespan (that is, the coronal section dataset), we used a correlation plot visualization consisting of a two-dimensional histogram of cell frequencies across bins defined by predicted age and actual age that is visualized as a heat map, a scatter plot of the median predicted age of cells across each sample as a function of the actual age, and a line of best fit for the median predicted ages as a function of actual age is shown in black and computed using numpy.polyfit with deg = 1. This type of visualization emphasizes the quality of median predicted age at the sample level across many different actual age values. Generally, in this visualization, the range of predicted ages will be larger than the range of actual ages due to heterogeneity in the predicted ages but not in the actual ages of cells from the same mouse. This may be especially pronounced for highly abundant cell types like excitatory neurons. Predicted ages obtained from cross-validation using a leave-mouse-out approach will exhibit some regression to the mean age (that is, cells from older mice predicted to be younger and cells from younger mice to be predicted to be older) due to different mean actual ages in each of the cross-validation training datasets.
For datasets with bimodal distribution of ages (that is, the sagittal section dataset) or with three or fewer distinct age groups (that is, all external datasets), we used a density plot visualization of predicted ages of cells across different age groups or different experimental conditions. A kernel density estimate was constructed for each group of predicted ages using seaborn.kdeplot with default settings. This type of visualization emphasizes the distribution of predicted ages in a small number of groups at the cell level. For the coronal sections dataset and sagittal sections dataset, we used both types of representations.
For comparisons of clock performances, we compare the mean absolute error and Pearson correlation between predicted age and actual age either with respect to the main spatial ageing clocks or between different data subsets or clocks to make our conclusions.
Application of spatial ageing clocks on external datasets
To apply the spatial ageing clocks to predict age for cell gene expression profiles in external spatial transcriptomics datasets, we applied the following general procedure. First, we filtered the external dataset to only include genes present in our MERFISH panel of 300 genes and only cell types also represented among the spatial ageing clocks. Then, we normalized and log-transformed the raw expression values using the same approach as for our MERFISH data and applied SpatialSmooth (\(\alpha =0.8\)) to the log-normalized values for each cell type independently. For clock genes that are not present in the external dataset, we use the training data for the clock as reference data in the SpaGE algorithm104 (n_pv = 15) to impute the expression of the missing genes. Negative imputed values were clipped to zero. We performed imputation for genes from our MERFISH panel that were missing from these datasets (228 genes in the 140-gene MERFISH coronal section dataset, 5 genes in the single-nuclei RNA-seq dataset, 36 genes in the single-cell RNA-seq dataset, 225 genes in the LPS dataset, 128 genes in the Alzheimer’s mouse model dataset, 240 genes in the global demyelination through the EAE dataset, 236 genes in the localized demyelination injury dataset). Finally, for each cell type, we applied the corresponding spatial ageing clock to generate predicted ages from the smoothed spatial gene expression values. For single-nuclei RNA-seq datasets, which lack spatial information, we used the pseudobulk approach from a previous model4 with 20 cells contributing to each pseudocell instead of SpatialSmooth.
For all evaluations, we also quantified the magnitude of difference in median predicted age (in units of months) and the 95% confidence interval for this difference, which was computed from the empirical distribution of differences in median predicted ages using 1,000 bootstrap samples of individual mice within each condition. These statistics are reported in Supplementary Table 12. In some cases, imputation resulted in biased age predictions, but the differences in predicted age across ages and conditions were generally robust. In applying spatial ageing clocks with imputation, we recommend comparing the predicted age to known ages in the dataset to calibrate interpretations. Generally, we observed lower spatial ageing clock performance for cell types with low transcriptomic changes with age such as neurons (see Fig. 1e) or those with limited marker genes in external datasets such as neuroblasts (see Fig. 2d). Results were also generally consistent for all applications of the spatial ageing clocks without SpaGE imputation and when using spatial ageing clocks trained with the 220 genes (see Supplementary Table 7 for list) with less than 5% estimated spillover based on an internal Vizgen metric using gene expression variations influenced by local cellular composition (Supplementary Table 12).
Age acceleration calculation
We computed age acceleration for the predicted age of each cell to measure the deviation from its expected predicted age (that is, the average of predicted age across all cells from a given cell type and mouse). For each sample \(k\in \{1,\ldots ,K\}\) and cell type \(p\in \{1,\ldots ,P\}\), we define the set of cells belonging to both sample \(k\) and cell type \(p\) as \(Q\), and the age acceleration for each cell \(i\in Q\) is defined as:
$${(\text{Age acceleration})}_{i}={(\text{Predicted age})}_{i}-E\,\{{(\text{Predicted age})}_{j}|\,j\in Q\}$$
Single-nuclei RNA-seq ageing data processing
We applied our spatial ageing clocks to publicly available single-nuclei RNA-seq data on the CTX and STR of juvenile (0.93 months) and old (20.93 months) female C57BL/6 J mice17. We downloaded processed data objects containing the scaled log-normalized gene expression from https://cellxgene.cziscience.com/collections/31937775-0602-4e52-a799-b6acdd2bac2e and mapped several cell types to our cell-type classifications. We modified the preprocessing and imputation steps outlined in ‘Application of spatial ageing clocks on external datasets’ to account for the scaled log-normalized expression being used as input. Predicted ages were highly consistent when using a scaled log-normalized coronal section dataset for imputation.
Single-cell RNA-seq ageing data processing
We applied our spatial ageing clocks to publicly available single-cell RNA-seq data on whole-brain tissue (without hindbrain regions) of young (2–3 months) and old (21–22 months) male C57BL/6 J mice5. We downloaded processed datasets from https://portals.broadinstitute.org/single_cell/study/aging-mouse-brain and mapped several cell types to our cell-type classifications.
140-gene MERFISH ageing data processing
We applied our spatial ageing clocks to our previously published 140-gene MERFISH spatial transcriptomics data on whole-brain coronal sections of young (0.93 months), middle-aged (5.58 months) and old (20.93 months) male C57BL/6 JN mice29. We mapped several cell types to our cell-type classifications.
MERFISH ageing and LPS data processing
We applied our spatial ageing clocks to publicly available MERFISH spatial transcriptomics data on the CTX and STR of female C57BL/6 J mice in juvenile (0.93 months), young (5.58 months), old (20.93 months), and lipopolysaccharide (LPS)-injected conditions17. We downloaded processed data objects containing the scaled log-normalized gene expression from https://cellxgene.cziscience.com/collections/31937775-0602-4e52-a799-b6acdd2bac2e and mapped several cell types to our cell-type classifications and used existing region annotations. We modified the preprocessing steps outlined in ‘Application of spatial ageing clocks on external datasets’ to account for the scaled log-normalized expression being used as input. Predicted ages were highly consistent when using a scaled log-normalized coronal section dataset for imputation. We included all control (ageing) and LPS-injected mice in our analyses.
STARmap Alzheimer’s disease data processing
We applied our spatial ageing clocks to publicly available STARmap PLUS spatial transcriptomics data on the CTX and hippocampus of male TauPS2APP (Alzheimer’s disease model) and non-transgenic control mice across two ages (8 and 13 months)22. We downloaded processed data from https://singlecell.broadinstitute.org/single_cell/study/SCP1375 and mapped several cell types to our cell-type classifications and used existing region annotations. We included all control and Alzheimer’s disease mice in our analyses. The percentage of zero counts was higher in this dataset (94%) than in the coronal section dataset used for training the spatial ageing clocks (70%).
In situ sequencing global demyelination data processing
We applied our spatial ageing clocks to publicly available in situ sequencing spatial transcriptomics data on whole-brain coronal sections of male and female C57BL/6J mice (2.5 months) across global demyelination and control conditions23. Global demyelination models consisted of induction of EAE via injection of myelin oligodendrocyte glycoprotein. We downloaded processed data from https://zenodo.org/records/8037425 and selected only whole-brain coronal sections. We mapped several cell types to our cell-type classifications and used existing region annotations. The percentage of zero counts in this dataset (76%) was similar to that in the coronal section dataset used for training the spatial ageing clocks (70%).
MERFISH localized demyelination injury data processing
We applied our spatial ageing clocks to publicly available MERFISH spatial transcriptomics data of three whole-brain coronal sections at different depths collected from young (3–4 months) male C57BL/6J mice that were subjected to demyelination injury via stereotactic injection of lysophosphatidylcholine at coordinates (from bregma): (X, ±1.0 mm; Y, −0.1 mm)21. We downloaded processed data and metadata from Gene Expression Omnibus (GSE202638) and mapped several cell types to our cell-type classifications. Due to the spatially localized nature of the demyelination injury and the lack of control conditions in this dataset, we analysed the predicted ages of cells in this dataset by spatially visualizing all cells by their positive age acceleration (negative values floored at zero) across each of the coronal sections to determine spatial patterns in age acceleration with respect to the site of injury for all cell types and for different major cell types. The percentage of zero counts was higher in this dataset (94%) than in the coronal section dataset used for training the spatial ageing clocks (70%).
Effect of interventions on predicted age
For a given cell type and experiment, to quantify the difference in predicted ages between two experimental conditions, we computed the difference in median predicted ages between cells of the two conditions. Specifically for interventions, we subtracted the median predicted age of cells belonging to mice in the control condition (sedentary, control) from the median predicted age of cells belonging to mice in the intervention condition (exercise, OSKM, LPS or Alzheimer’s disease). Positive values indicate that the intervention has an accelerated ageing effect and negative values indicate that the intervention has a rejuvenating effect on the cells of the mice. For the whole-body partial reprogramming experiment, since there was a difference in the mean age of mice in the control and OSKM conditions of 0.3 months, we corrected for this difference by adding an intercept to the predicted ages before computing the effect. We computed this effect at the global levels (all cells of a cell type) and at the regional level (all cells of a cell type within a defined anatomic region). For the adverse interventions datasets (LPS and Alzheimer’s disease), we used the existing anatomic region annotations and mapped them to the closest region label in our study. For all global comparisons, we quantified the difference in median predicted age (in units of months) and the 95% confidence interval for this difference, which was computed from the empirical distribution of differences in median predicted ages using 1,000 bootstrap samples of individual mice within each condition. These statistics are reported in Supplementary Table 12.
Spatial visualizations of the effect of interventions on predicted age involved computing an effect value for each cell in a sample in the intervention condition. Specifically, for each cell in the intervention condition, we subtracted the median predicted age of cells belonging to mice in the age-matched control condition (sedentary, control) from the predicted age of that cell, and then visualized that cell by its spatial coordinates and coloured by the computed value.
Cell proximity effects of ageing and rejuvenation
The cell proximity effect measures the effect that an effector cell type has on the transcriptomic ageing of a target cell type by comparing target cells near effector cells to those that are far from effector cells. To compute the distribution of nearest neighbour distances, we constructed a triangulation mesh graph connecting neighbouring cells on a given sample using squidpy.gr.spatial_neighbors with delaunay=True103. We used the centroid of each cell to compute distances between cells. We computed subregion-specific distance cutoffs for calling nearby cells as the average of the median neighbour-neighbour distances measured across all samples. The distance cutoffs (in micrometres) for the ageing (coronal) study were: CC/ACO: 24.89, CTX_L1/MEN: 25.91, CTX_L2/3: 24.05, CTX_L4/5/6: 27.24, STR_CP/ACB: 21.65, STR_LS/NDB: 20.36, VEN: 17.86. The distance cutoffs (in micrometres) for the exercise study were: CC/ACO: 23.58, CTX_L1/MEN: 22.13, CTX_L2/3: 21.80, CTX_L4/5/6: 24.81, STR_CP/ACB: 20.75, STR_LS/NDB: 19.82, VEN: 16.23. Using these subregion-specific distance cutoffs, we identified target cells near effector cells (‘near’) and matched them to target cells far from any effector cells (‘far’) for a given target cell type and effector cell type. First, we computed for each cell, the shortest Euclidean distance to any effector cell in the same sample for each effector cell type. Then, for each sample and combination of target cell type and effector cell type, we labelled all target cells with shortest Euclidean distance to effector cell type that was less than the corresponding distance cutoff as ‘near’ and match them with ‘far’ target cells in the same subregion and sample that were furthest away from the effector cell type with shortest distance greater than the cutoff for ‘near’ cells. After obtaining matched sets of ‘near’ and ‘far’ target cells across all samples, we combined them into a single set to estimate the proximity effect of the effector cell type on the target cell type. The proximity effect was computed as Cohen’s d measure of effect size between the age acceleration of the ‘near’ target cells compared with the age acceleration of the ‘far’ target cells. The variance in age acceleration was similar between ‘near’ and ‘far’ target cell groups (Extended Data Fig. 9a). We filtered out comparisons with less than 50 ‘near’ cells or less than 50 ‘far’ cells. Positive proximity effect values indicate a pro-ageing effect exerted by the effector cell type on the target cell type and negative proximity effect values indicate a pro-rejuvenating effect exerted by the effector cell type on the target cell type. We also computed the normalized frequency of proximity interactions as the number of ‘near’ target cells divided by the total number of target cells in the study and the statistics and P value from the associated two-sided Student’s t-test between the age acceleration of the ‘near’ and ‘far’ target cells.
We implemented several variations to the standard proximity effect analysis. For region-specific proximity effects, we restricted the proximity effect analysis to only cells within each of the anatomic regions. For cell proximity effects using alternative definitions of ‘far’ cells, we used two orthogonal approaches. The first approach consisted of random sampling of ‘far’ cells from all target cells with shortest distance greater than the cutoff for ‘near’ cells in the same subregion and sample. The second approach consisted of selecting the set of ‘far’ cells with total raw transcript counts closest to the mean total raw transcript count of the ‘near’ cells in the same subregion and sample. For cell proximity effects based on predicted ages obtained from non-spatial ageing clocks such as spatial ageing clocks that did not use the SpatialSmooth step for prediction (SingleCell (SS)) or single-cell ageing clocks using cell-type-specific pseudobulking of gene expression (1,000 bootstrap samples of 30 cells following prescribed procedures4) during training and prediction (SingleCell (PB)), we used predicted ages from these aforementioned clocks in lieu of predicted ages from the spatial ageing clocks and performed cell proximity analysis following the original setup. For the ‘spillover filtered’ cell proximity effects, we used predicted ages from the spatial ageing clocks trained on the 220 genes with less than 5% estimated spillover based on an internal Vizgen metric using gene expression variations influenced by local cellular composition. For cell proximity effects on datasets other than the coronal section dataset, we computed dataset-specific distance cutoffs for subregions and regions and computed proximity effects using the original setup for all cell types that were mapped to our cell-type classifications. We calculated separate proximity effects for cells across the entire external dataset and for cells restricted to the control conditions. To assess the effect of interventions on cell proximity effects, we computed cell proximity effects separately for each experimental condition. We excluded the 140-gene MERFISH dataset from proximity effect analysis because this dataset was spatially sparse with distance cutoffs that were several times greater than that of other datasets, probably owing to lower cell detection/segmentation frequency.
The average proximity effect for each effector cell type was computed as the average proximity effect across all target cell types after filtering. To assess the most impacted target cell types, the mean absolute proximity effect for each target cell type was computed by averaging the absolute value of the proximity effect across all effector cell types after filtering.
Spatial permutation analysis
Spatial permutations can be used to assess frequency of cell–cell interactions occurring by chance due to regional differences in cell-type proportion105. We performed subregional spatial permutations of cells as a negative control for the spatial proximity effects measured for effector cells. In each permutation, for each sample in the ageing (coronal) study, we randomly permuted the spatial coordinates of each cell using numpy.random.permutation. This permutation was performed separately for each subregion to maintain regional heterogeneity. Then, we measured proximity effects on the permuted dataset. We repeated these permutations twenty times using random seeds drawn uniformly from 0 to 50,000 with a generating random seed of 444. We computed the median and 95% confidence interval for the proximity effect of each effector cell type across these twenty permutations.
Spatial visualization of age acceleration near effector cells
We generated visualizations of cells coloured by their age acceleration around key effector cell types (for example, T cells and NSCs) by selecting an effector cell, drawing a square bounding box centred on the effector cell with edge lengths of 200 μm, and then visualizing all cells with centre coordinates inside of the bounding box. We manually selected visualizations from a middle-aged (18.8 months) mouse for T cells and NSCs.
Impact of potential spillover transcripts on proximity effect
We verified that the cell proximity effects were unlikely to be confounded by spillover of transcripts between nearby effector and target cells arising from segmentation errors, which can be particularly prominent for cells near NSCs in the densely packed VEN region (Supplementary Fig. 7a). In a first approach to test this, we computed proximity effects based on predicted ages obtained from spatial ageing clocks trained on a subset of genes after filtering out genes with high estimated spillover rate (see Methods) (Supplementary Fig. 7b). In a second approach to test this, we performed area-restricted proximity effect analysis, which excluded cells within a small radius of the effector cells (see ‘Area-restricted cell proximity effects’), for spatial ageing clocks and two non-spatial ageing clocks (Supplementary Fig. 8a,b). We verified that this area-restricted approach drastically attenuated transcript spillover effects in all target cell types near T cells and NSCs (Supplementary Fig. 8c). Both approaches corroborate that T cells have the most pro-ageing proximity effect and NSCs have the most pro-rejuvenating proximity effect (Supplementary Figs. 7b and 8b).
Area-restricted cell proximity effects
To remove the potential influence of transcript spillover from cell segmentations on the cell proximity effect analysis, we developed an alternative approach to compute area-restricted cell proximity effects, where ‘near’ cells are labelled using two cutoff distances instead of the single cutoff distance in the standard proximity effect analysis (see Supplementary Fig. 8a). We set the larger of the two cutoff distances equal to double the subregion-specific distance cutoff from ‘Cell proximity effects of ageing and rejuvenation’. We set the smaller cutoff distance to be 15 μm less than the larger cutoff distance. We label target cells with shortest Euclidean distance to effector cell type that is greater than the smaller distance cutoff but less than the larger distance cutoff as ‘near’ and match them with ‘far’ target cells in the same region and sample that are farthest away from the effector cell type with shortest distance greater than the larger distance cutoff for ‘near’ cells. Then, area-restricted cell proximity effects are calculated using the same approach as described in ‘Cell proximity effects of ageing and rejuvenation’ except with these modified ‘near’ and ‘far’ cell labels.
Spatial gradation of cell proximity effects
We computed cell proximity effects as a function of the unit distance between the effector and target cells, where the unit distance is defined as a scalar value that is multiplied by the subregion-specific distance cutoffs from ‘Cell proximity effects of ageing and rejuvenation’. These unit distances are used as new cutoff distances to compute cell proximity effects. We also analysed spatial gradation of area-restricted cell proximity effects by setting the larger of the two cutoff distances to the unit distance and the smaller cutoff distance to be 15 μm less than the larger cutoff distance, which ensured no overlap of cells defined as ‘near’ between integer multiples of the cutoff distances.
Activation/inflammation glia signatures and subtype identification
We computed microglia activation scores and oligodendrocyte/OPC inflammation scores by summing the log-normalized expression of all shared genes between gene signature sets curated from several published sets17,106,107 and the 300 genes in the MERFISH dataset. The microglia activation score consisted of gene expression from Apoe, B2m, C1qa, Cd69, Cd9, Il1a, Il1b, Il6 and Lyz2. The oligodendrocyte/OPC inflammation score consisted of gene expression from C4b, Cdkn1a, H2-D1, Ifit1 and Stat1.
To control for activation/inflammation status in the T cell proximity effect, we selected the cutoff for classifying microglia (activated) and oligodendrocyte (inflamed) subtypes as the top 0.2% highest scores since that was the highest percentage cutoff for which there was no statistically significant difference in the distribution of activation/inflammation scores between activated/inflamed cells that were near or far from T cells. We performed all associated proximity effect analysis using the normal and activated/inflamed subtype labels for microglia and oligodendrocytes.
Cell-type perturbation modelling with deep learning
Deep learning approaches such as GNNs can be leveraged to predict the effects of in silico perturbations. For each section in the coronal section dataset, we constructed a global graph connecting neighbouring cells on a given sample using squidpy.gr.spatial_neighbors with delaunay=True103 and pruned edges connecting neighbouring cells with distance greater than 200 μm. To define local cell graphs, we extracted induced 2-hop subgraphs of the global graph by random sampling of at most 100 centre cells per cell type for T cells and NSCs. To increase the heterogeneity of these local cell graphs, we further augmented these local cell graphs by inducing two-hop subgraphs centred on all cells within the first set of local subgraphs. For each local cell graph, we defined the neighbourhood ageing as the average age acceleration of all cells in the graph. We defined cell (node) features as a one-hot vector for cell type using the 18 cell-type annotations. We trained a GNN model using PyTorch Geometric108 to predict neighbourhood ageing from the features of a local cell graph. The GNN model consisted of a three-layer graph isomorphism network (GIN), with node embedding updates for each layer modelled by a linear transformation with hidden dimensionality of 16 followed by batch normalization, ReLU transformation, and a final linear transformation with hidden dimensionality of 16. The first two GIN layers were followed by a ReLU transformation and the output of the final GIN layer was globally pooled by addition before linearly transformed to predict neighbourhood ageing. We trained the GNN model using a balanced mean-squared error loss109, the Adam optimizer110 with a learning rate of 0.0001, and a batch size of 64 for 50 epochs.
To model ‘loss-of-function’ perturbations for a given local cell graph, we mutated the centre cell of the graph into a random cell type drawn from a uniform distribution across all cell types in the dataset excluding the original cell type. We performed ‘loss-of-function’ perturbations for local cell graphs that had T cells and NSCs as centre cells. To model ‘gain-of-function’ perturbations for a given local cell graph, we replaced the centre cell of the graph with the specified cell type by modifying its node features. We performed ‘gain-of-function’ perturbations for all local cell graphs and mutated the centre cells to either T cells or NSCs. We excluded ‘gain-of-function’ perturbations resulting in unperturbed local cell graphs from our analysis (for example, graph with T cell as centre cell replaced by T cell). For both ‘loss of function’ and ‘gain-of-function’ experiments, we used endothelial cells as a negative control setting. Using the GNN model, we evaluated the effect of both ‘loss-of-function’ and ‘gain-of-function’ perturbations by predicting the neighbourhood ageing for the unperturbed local cell graph and the neighbourhood ageing for the perturbed local cell graph. The effect on neighbourhood ageing was then measured as the unperturbed neighbourhood ageing subtracted from the perturbed neighbourhood ageing (positive values indicating a pro-ageing perturbation and negative values indicating a pro-rejuvenating perturbation).
TISSUE imputation
To augment the 300 genes measured through MERFISH to more than 12,000 genes, we performed uncertainty-aware spatial gene expression imputation, using a method we recently developed, TISSUE29, as a wrapper around the SpaGE104 and Tangram111 imputation algorithms. TISSUE is an algorithm that provides an uncertainty-aware framework for performing differential gene expression and signature analysis on imputed spatial gene expression with marked reductions in false discovery rates compared with other approaches29. This imputation approach consists of jointly mapping transcriptomes from our MERFISH study with single-cell RNA-seq datasets collected from mice across multiple ages and containing the brain regions of interest (for example, VEN, CC or STR)3,4, followed by prediction of new gene expression profiles using the single-cell RNA-seq dataset as reference and differential analysis of gene expression and gene signatures (Extended Data Fig. 12a). The log-normalized gene expression data from the first dissociated single-cell RNA-seq dataset3 were used for imputation of gene expression for T cells (8,170 genes imputed in total), and log-normalized gene expression data from the second dissociated single-cell RNA-seq dataset4 were used for imputation of gene expression for all other cell types (12,719 genes imputed in total). We applied the TISSUE algorithm on a subset of the coronal section dataset (see below for details). We used the raw MERFISH counts normalized by cell segmentation volume as input for imputation. We used TISSUE with the SpaGE imputation algorithm104 and Tangram imputation algorithm111, under default settings such as 10 folds of cross-validation, 4 stratified gene groups, 1 stratified cell group, and 100 multiple imputations (for hypothesis testing purposes). Evaluation of TISSUE calibration quality and SpaGE/Tangram imputation performance were conducted using the TISSUE software package and associated code29. Using cross-validation, we verified that the imputed gene expression values were positively correlated with the actual gene expression values and generally had small absolute prediction errors (Extended Data Fig. 12b).
We filtered cells in the coronal section dataset to include NSCs and T cells, and all target cells labelled ‘near’ or ‘far’ with respect to NSCs or T cells under the proximity effect framework. We re-centred each sample to grid lattice points such that no samples spatially overlap. In addition to ‘near’ and ‘far’ labels, we also labelled each T cell and NSC by the strength of their proximity effects. First, we computed the neighbourhood age acceleration as the average age acceleration of all cells within the maximum subregion distance cutoff for each T cell or NSC. Then, we labelled all T cells with the 50% highest neighbourhood age acceleration as ‘more pro-ageing’ and the remaining T cells as ‘less pro-ageing’, and we labelled all NSCs with the 50% lowest neighbourhood age acceleration as ‘more pro-rejuvenating’ and the remaining NSCs as ‘less pro-rejuvenating’.
TISSUE differential gene expression analysis and pathway enrichment
To uncover potential mediators of T cell and NSC proximity effects, we used a two-step approach. First, we performed unbiased differential gene expression analysis using TISSUE29 on more than 12,000 imputed genes, with the goal of identifying common genes that were upregulated in cells near key effector cell types (T cells or NSCs) compared with cells far from the same effector cell type and matched by cell type and subregion (that is, the same procedure as for calculating proximity effects). We verified that cells near T cells or NSCs did not include a cell type that was substantially more represented than all other cell types and that all 14 cell types with spatial ageing clocks were present (Supplementary Table 16), and we checked that the variances in imputed gene expression were similar between the comparison groups (Extended Data Fig. 12c). In the second step, to better understand potential mechanisms of action, we performed targeted statistical comparisons of imputed gene signatures using TISSUE and gene expression in the MERFISH datasets for processes related to the most enriched terms from the unbiased analyses in both effector cell types (T cells or NSCs) and target cells (nearby cells).
We used TISSUE to perform uncertainty-aware differential gene expression analysis on the imputed gene expression. TISSUE relies on a multiple imputation framework with statistical guarantees that are only extensible to t-tests (non-parametric TISSUE tests rely on a less rigorously defined P value transformation)29. Given the variability in imputation quality, we chose the TISSUE framework, which was well-suited for this task despite some limitations of the t-test. We performed two sets of comparisons using TISSUE two-sided t-test for hypothesis testing of gene signatures. First, we compared gene signature scores between ‘near’ and ‘far’ target cells with respect to either T cells or NSCs. Second, we compared gene signature scores between more pro-ageing and less pro-ageing T cells or between more pro-rejuvenating and less pro-rejuvenating NSCs. We determined significant DEGs as genes with permissive cutoff of P < 0.05 from TISSUE two-sample t-tests for both SpaGE and Tangram imputed expression to reduce false discoveries resulting from variability or biases in imputation quality. We further filtered these genes such that the associated TISSUE t-statistic were of the same sign across the tests on SpaGE and Tangram imputations. We did not perform any filtering of genes by the log fold change in imputed gene expression across conditions because unlike the TISSUE P value, the log fold change is not calibrated for uncertainty in spatial gene imputation.
GO biological process (2023) gene set enrichment was performed using the EnrichR framework112 accessed through gseapy113 (version 1.0.4). We separated DEGs for each comparison into upregulated DEGs (positive TISSUE t-statistic) and downregulated DEGs (negative TISSUE t-statistic) and used these gene lists as inputs into EnrichR to obtain enrichment statistics for different gene sets. We selected the top five gene sets for each comparison, which were generally representative of all significantly enriched gene sets.
Test for equal variances between groups
We performed Levene’s test for equal variances using the scipy.stats.levene implementation with default settings (center = median, proportiontocut = 0.05) for analyses relying on Cohen’s d or variations of the Student’s t-test. Levene’s test was applied to test for equal variance of age acceleration of ‘near’ and ‘far’ cell groupings for each cell proximity effect. Levene’s test was also applied to test for equal variance of TISSUE (SpaGE) imputed gene expression of ‘near’ and ‘far’ cell groupings with respect to T cells and NSCs. In both cases, the number of samples between the compared groups was equal, a setting in which Cohen’s d and Student’s t-test are usually robust against unequal variances across groups.
TISSUE gene signature scores
We further modified the hypothesis testing framework in TISSUE to perform testing for gene signature scores (sum of imputed gene expression values across all genes in a GO gene set). We have made this modification publicly available in the TISSUE package29 as tissue.downstream.multiple_imputation_gene_signature and we set n_imputations=100 for all gene signature tests. For imputed gene signature scores, we performed the same two sets of comparisons as outlined in ‘TISSUE differential gene expression analysis and pathway enrichment’. For similar comparisons involving expression of individual genes (for example, Bst2, Stat1, Cd9 and Vegfa) instead of gene signatures, we used the measured (log-normalized) gene expression values in the MERFISH dataset instead of imputed values.
Immunofluorescence staining of brain sections
All immunostainings were performed on male C57BL/6JN mice at the indicated ages. Mice were sedated with 0.8 ml 2.5% vol/vol Avertin (Sigma-Aldrich, T48402-25G) in PBS (Corning, 21-040-CV) and perfused via the left ventricle of the heart with 5 ml of PBS (Corning, 21-040-CV) with heparin sodium salt (50 U ml−1, Sigma-Aldrich, H3149-50KU) to remove circulating blood cells followed by 4% paraformaldehyde (Electron Microscopy Sciences, 15714) in PBS. Brains were post-fixed overnight in 4% paraformaldehyde (Electron Microscopy Sciences, 15714) and then dehydrated in 30% sucrose (Sigma-Aldrich, S3929) for 72 h. Brains were embedded in OCT compound (Fisher Healthcare Tissue Plus, 4585), sectioned in 20-μm coronal sections using a cryostat (Leica, CM3050S), and then mounted on electrostatic glass slides (Fisher Scientific, 12-550-15). Coronal sections at a similar depth to the MERFISH coronal sections were collected. For the following steps, all brain sections were processed simultaneously within each experiment. For immunofluorescence staining, brain sections were washed with PBS for 5 min, permeabilized in ice-cold methanol with 0.1% Triton X-100 (Fisher Scientific, BP151) for 15 min, and washed 3 times with PBS for 5 min. We performed antigen retrieval by placing the brain sections in 10 mM sodium citrate buffer (pH 6.0; 2.94 g Tri-sodium citrate dihydrate (Sigma-Aldrich, S1804) in 1,000 ml milliQ H2O adjusted to pH 6.0 with 1 N HCl) + 0.05% Tween-20 (Sigma-Aldrich, P1379-1L) at 85 °C in a water bath for 2 h. Brain sections in buffer were cooled to room temperature for 20 min and then washed twice with PBS for 3 min. Sections were blocked for 30 min at room temperature with block buffer consisting of 5% normal donkey serum (ImmunoReagents, SP-072-VX10) and 1% BSA (Sigma, A7979) in PBS. Primary antibodies were diluted in blocking buffer and primary antibody staining was performed overnight at 4 °C. Primary antibodies used were as follows: anti-GFAP (1:1,000 dilution, Abcam, ab53554), anti-Ki67 (1:500 dilution, Thermo Fisher Scientific, 14-5698-82), anti-DCX (1:500 dilution, Millipore Sigma, AB2253), anti-EGFR (1:200 dilution, Millipore Sigma, 06-847), anti-STAT1 (1:500 dilution, Cell Signaling Technology, 14994), anti-CD3 (1:500, Abcam, ab11089), anti-S100A6 (1:500, Abcam, ab181975), anti-CD9 (1:100, Thermo Fisher, eBioKMC8 (KMC8)), anti-CPT1A (1:200, Abcam, ab128568). Sections were washed 3 times with PBS and 0.2% Tween-20 for 10 min at room temperature followed by twice with PBS for 15 min. Secondary antibodies were diluted in blocking buffer and secondary antibody staining was performed for 2 h at room temperature. Secondary antibodies used were as follows: Donkey anti-Goat 647 (1:500 dilution, Invitrogen, A21447), Donkey anti-Rabbit 488 (1:500 dilution, Invitrogen, A-21206), Donkey anti-Guinea pig 594 (1:500 dilution, Jackson ImmunoResearch, 706-585-148), Donkey anti-Rat 647 (1:500 dilution, Invitrogen, A48272), Donkey anti-Rat 488 (1:500, Invitrogen, A21208), Donkey anti-Mouse 647 (1:500, Invitrogen, A31571), Donkey anti-Rabbit 568 (1:500, Invitrogen, A10042). DAPI (1:500, Thermo Fisher, 62248) was added during secondary antibody staining. Sections were washed 3 times with PBS and 0.2% Tween-20 for 10 min followed by 3 times with PBS for 5 min. Sections were mounted with ProLong Gold Antifade Mountant with DAPI (Thermo Fisher, P36931) and a coverslip.
For immunofluorescence experiments corresponding to spatiotemporal marker expression, images were acquired on a Zeiss LSM 900 with Zeiss ZEN Blue 3.0 software or a Zeiss LSM 980 confocal microscope with Zeiss ZEN Blue 3.3 software using a 10× objective and automatic tiling of entire coronal brain sections. The same equipment and microscope acquisition settings were used for different brain sections with the same antibody panel. Tile images were stitched using Zeiss ZEN Blue software. Image brightness and contrast were adjusted in ImageJ (1.53n) to enhance visualization with the same settings applied to all images shown for each antibody staining panel.
For immunofluorescence experiments corresponding to T cell and NSC proximity effects, images were acquired with a 60× objective on a Nikon Eclipse Ti confocal microscope equipped with a Zyla sCMOS camera (Andor) and NIS-Elements software (AR 4.30.02, 64-bit). For immunofluorescence experiment corresponding to T cell proximity effects, we acquired at least two images containing T cells and at least one image without any T cells from each of four anatomic regions (CC, CTX, striatum and adjacent regions (STR) and VEN) for each brain section. The same image acquisition settings were used for all brain sections in this experiment. For immunofluorescence experiment corresponding to NSC proximity effects, we acquired individual images tiling the entire right and left VEN for each brain section. The same image acquisition settings were used for all brain sections in this experiment.
Immunofluorescence quantification of T cell proximity effect via interferon response
Immunofluorescence imaging was performed on brain sections from old (28 months) mice. All cell segmentations and cell-type annotations were performed using automated pipelines in QuPath 0.5.1. For all images, cell nuclei were automatically segmented based on DAPI intensity, and then nuclear segmentation masks were expanded by 2 μm to define the cell segmentations. Cells were labelled as CD3+ using a manually determined threshold for mean cell CD3 intensity. The same threshold was used across all images and CD3+ cells were annotated as T cells. For all CD3– cells, we defined cells as ‘near’ if they were located within 50 μm of a T cell based on the centroids of cell segmentations and otherwise defined them as ‘not near’. For each anatomic region (CC, CTX, STR and VEN) and each cell proximity definition (‘near’ or ‘not near’), we quantified the mean STAT1 intensity by averaging the mean cell STAT1 intensity across each section and then by averaging these values across each mouse (3–5 coronal sections per mouse, 8 mice per condition). We were unable to discover any T cells in the STR region for one of the eight mice due to low CD3 intensity in its brain sections. The mean STAT1 intensity values were normalized for each independent experiment by dividing them by the corresponding mean ‘not near’ STAT1 intensity.
Immunofluorescence quantification of NSC proximity effect via exosomes and fatty acid oxidation pathway
Immunofluorescence imaging was performed on brain sections from young (3.5 months) mice. All cell segmentations and cell-type annotations were performed using automated pipelines in QuPath 0.5.1. For all images, regions of interest were defined along the lining of the VEN. Within these regions of interest, cell nuclei were segmented based on DAPI intensity, and then nuclear segmentation masks were expanded by 2 μm to define the cell segmentations. Cells were labelled as S100A6+ using a manually determined threshold for mean nuclear S100A6 intensity. The same threshold was used across all images and S100A6+ cells were annotated as NSCs and other cells were annotated as non-NSCs. To compare CD9 intensity between S100A6+ and S100A6− cells, we quantified the mean CD9 intensity by averaging the mean cell CD9 intensity across each section and then by averaging these values across each mouse (three coronal sections per mouse, five mice per condition). The mean CD9 intensity values were normalized for each independent experiment by dividing them by the CD9 intensity in S100A6− cells. For all cells, we defined cells as ‘near’ if they were located within 20 μm of a NSC (excluding itself) based on the centroids of cell segmentations and otherwise defined them as ‘not near’. To compare CPT1A intensity between ‘near’ and ‘not near’ cells, we quantified the mean CPT1A intensity by averaging the mean cell CPT1A intensity across each section and then by averaging these values across each mouse (three coronal sections per mouse, five mice per condition). The mean CPT1A intensity values were normalized for each independent experiment by dividing them by the corresponding mean ‘not near’ CPT1A intensity. To examine the correlation between CPT1A intensity of all ‘near’ cells and the CD9 intensity of the nearest NSC, we matched all ‘near’ cells to their nearest NSC defined by the minimum non-zero Euclidean distance between centroids and computed the Pearson and Spearman correlation between the paired CPT1A and CD9 intensities. To check the robustness of the correlation analysis to technical variability, we also re-imaged all sections using a Zeiss LSM 900 with Zeiss ZEN Blue 3.0 software and performed the same image analysis. Results were highly consistent with the reported findings.
Statistics and reproducibility
For measures of Pearson and Spearman correlation and their associated test statistics (P values and 95% confidence intervals), we used the implementation from scipy.stats.pearsonr and scipy.stats.spearmanr. For most comparisons, we used the two-sided Mann–Whitney tests according to the implementation from scipy.stats.mannwhitneyu. We used Cohen’s d for cell proximity effect analysis in which similar variances in age acceleration were typically observed between ‘near’ and ‘far’ cell groups (Extended Data Fig. 9a). GO biological process enrichment analysis was performed with topGO101 (R package version 2.54.0) using Fisher’s exact test for all Biological Process terms and using all other genes measured by MERFISH as background. For other data analysis and plotting tasks, we used Python (3.8.13), pandas (1.4.2), matplotlib (3.5.1), seaborn (0.12.2), numpy (1.21.6), scipy (1.8.0), sklearn (1.0.2), anndata (0.8.0), scanpy (1.9.1), squidpy (1.1.2), tissue-sc (0.0.2), tangram-sc (1.0.3), spage (accessed September 1, 2022), gseapy (1.0.4), umap-learn (0.5.3) and statsmodels (0.13.2).
For the gene signature and differential gene expression analysis with imputed spatial gene expression, we used the two-sample two-sided TISSUE t-test with specifications outlined in the ‘TISSUE imputation’ section. The TISSUE t-test is the only statistical test that considers uncertainty in imputed gene expression and to further decrease false discoveries, we selected genes with consistent differential expression across two imputation methods (SpaGE104 and Tangram111). For most comparisons using the TISSUE t-test, variances in imputed gene expression were similar across groups (Extended Data Fig. 12c). We did not perform multiple hypothesis correction of P values for the targeted gene signature analysis using TISSUE imputation because gene signatures were manually selected. Reported P values are from the SpaGE imputed gene expression. For pathway enrichment for DEGs identified with TISSUE and imputed gene expression, we used gseapy.enrichr with the ‘GO_Biological_Process_2023’ gene set, ‘mouse’ as the organism, and all imputed genes as background.
Technical reproducibility of MERFISH measurements was confirmed by profiling consecutive sections in the 140-gene MERFISH dataset29 (Extended Data Fig. 1a). For the MERFISH (300 gene) ageing dataset, we combined two cohorts (n = 10 mice for each cohort) and confirmed reproducibility of spatial ageing clock performance and generalization of key findings within each cohort (Extended Data Fig. 5b,c and Supplementary Fig. 6d,e). All experimental validation was performed on independent mice. Technical reproducibility of immunofluorescence image analysis for the NSC proximity effect mediation through exosomes and fatty acid oxidation was assessed with re-imaging of the same sections by a different experimentalist with a different microscope. Two separate immunofluorescence staining and imaging experiments for T cell proximity effect mediated through interferon response were performed and combined for analysis, with intensities normalized against each experiment independently.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
