You are viewing the site in preview mode

Skip to main content

Regional and aging-specific cellular architecture of non-human primate brains

Abstract

Background

Deciphering the functionality and dynamics of brain networks across different regions and age groups in non-human primates (NHPs) is crucial for understanding the evolution of human cognition as well as the processes underlying brain pathogenesis. However, systemic delineation of the cellular composition and molecular connections among multiple brain regions and their alterations induced by aging in NHPs remain largely unresolved.

Methods

In this study, we performed single-nucleus RNA sequencing on 39 samples collected from 10 brain regions of two young and two aged rhesus macaques using the DNBelab C4 system. Validation of protein expression of signatures specific to particular cell types, brain regions, and aging was conducted through a series of immunofluorescence and immunohistochemistry staining experiments. Loss-of-function experiments mediated by short hairpin RNA (shRNA) targeting two age-related genes (i.e., VSNL1 and HPCAL4) were performed in U251 glioma cells to verify their aging effects. Senescence-associated beta-galactosidase (SA-β-gal) staining and quantitative PCR (qPCR) of senescence marker genes were employed to assess cellular senescence in U251 cells.

Results

We have established a large-scale cell atlas encompassing over 330,000 cells for the rhesus macaque brain. Our analysis identified numerous gene expression signatures that were specific to particular cell types, subtypes, brain regions, and aging. These datasets greatly expand our knowledge of primate brain organization and highlight the potential involvement of specific molecular and cellular components in both the regionalization and functional integrity of the brain. Our analysis also disclosed extensive transcriptional alterations and cell–cell connections across brain regions in the aging macaques. Finally, by examining the heritability enrichment of human complex traits and diseases, we determined that neurological traits were significantly enriched in neuronal cells and multiple regions with aging-relevant gene expression signatures, while immune-related traits exhibited pronounced enrichment in microglia.

Conclusions

Taken together, our study presents a valuable resource for investigating the cellular and molecular architecture of the primate nervous system, thereby expanding our understanding of the mechanisms underlying brain function, aging, and disease.

Background

Owing to their genetic, physiological, and anatomical similarities to humans, non-human primates (NHPs) such as macaque, marmoset, sooty mangabey, African green monkey, and baboon have contributed enormously to human health research as models for investigating human diseases and assessing both the efficacy and safety of newly developed drugs [1,2,3]. They have been particularly widely used to investigate specific structures and functions in the brain [4,5,6], because of their complex brains and enhanced cognitive abilities compared to rodent models. The brain represents the most complex organ; it contains multiple distinct cell types within diverse regions and acts as a control center for behavior, emotions, and cognition [7, 8]. The precise modulation of neural circuitry among cell types and different brain regions is essential for robust brain functioning [9], while cellular and circuit dysfunction has the potential to cause brain disorders including autism, schizophrenia, and those conditions that often develop with advancing age such as Alzheimer’s disease (AD) and Parkinson’s disease (PD) [9,10,11,12,13,14].

Understanding brain function and brain disease requires the deep systematic characterization of different cells and their interactions across distinct brain regions. Rapid advances in single-cell RNA sequencing (scRNA-seq) technology have expanded our knowledge of the cellular composition of diverse tissues and organs [15,16,17,18], and multiple cell types have been found to characterize the structural organization of the vertebrate brain [5, 19,20,21,22,23]. However, collecting fresh tissue samples from many different brain regions is logistically challenging in primates, and particularly so in humans. In addition, the brain is difficult to process into single-cell suspensions and the morphology of brain cells is too irregular to withstand traditional single-cell approaches [24]. As a consequence, our knowledge and understanding of the cellular composition of different primate brain regions remains quite limited, despite an increasing number of scRNA-seq studies of primate brains from specific regions [5, 20,21,22,23]. Single-nucleus RNA sequencing (snRNA-seq), which can be applied to archived and frozen brain tissues, offers a compelling alternative. In this study, we have characterized the brain cellular transcriptome profiles of the rhesus macaque (Macaca mulatta) based on approximately 330,000 cells obtained from ten brain regions of young and aged individuals using snRNA-seq. Our work reveals the intricacy of the molecular and cellular organization of the primate brain and its dynamics with aging, providing an important resource for the in-depth analysis of the mechanisms underpinning brain function, aging and disease.

Methods

Animal

Rhesus macaques with no previously reported neuropsychiatric disorders were purchased from the Kunming Primate Research Center, Chinese Academy of Sciences (AAALAC accredited). The macaques were housed individually in cages and maintained under a 12-h light/dark cycle (from 7:00 AM to 7:00 PM), with humidity at 60% and a room temperature of 21 ± 2 °C. They received two daily servings of monkey chow and one serving of fruit, while water was available at all times. Their health was closely monitored by skilled veterinarians. The monkeys were anesthetized with ketamine (10 mg/kg, i.m.), followed by euthanasia with sodium pentobarbital (100 mg/kg, i.m.) and trans-cardiac perfusion with PBS. All experimental protocols had been approved by the Institutional Animal Care and Use Committee of the Kunming Institute of Zoology, Chinese Academy of Sciences (permit nos. IACUC-PE-2022-–01-003).

Sample collection for snRNA-seq

A total of 40 samples were collected from 10 brain regions of two young and two aged rhesus macaques. These specimens were obtained post-euthanasia, remaining from our previous study [25], and stored at − 80 °C in the Kunming Primate Research Center. The two young macaques included SY1 (female, 6 years old) and SY2 (male, 6 years old), while the two old macaques included SO1 (female, 17 years old) and SO2 (female, 24 years old). As we previously described [25], the brain regions were collected by a skilled technician according to a widely used macaque brain atlas (http://www.brainmaps.org) [26], including amygdala (AMY), putamen (PU), hippocampus (HIP), thalamus (TH), dorsolateral prefrontal cortex (DLPFC), cingulate gyrus (CG), superior temporal gyrus (STG), superior parietal lobule (SPL), visual cortex V4 (V4), and cerebellar cortex (CBC). Loci and histological sections of the brain regions were displayed in Additional file 1: Fig. S1. Before sample collection, all surgical instruments were sterilized prior to use. To prevent cross-contamination, fresh sterile scissors and tweezers were used for each sample. After dissection, the tissue samples were rinsed with RNAlater solution (AM7021, Ambion, USA) and then transferred into cryovials for immediate storage in liquid nitrogen. Bulk RNA-seq was conducted using parts of these samples in our previous study to investigate potential molecular mechanisms of brain aging in no-human primates [25]. For better understanding, we renamed the samples here and presented their old corresponding names used by our previous study [25] in Additional file 2: Table S1.

Nucleus isolation

snRNA-seq data were then obtained using the DNBelab C4 droplet-based platform by BGI-Shenzhen (Shenzhen, 518,083, China) as previously described [15, 27]. In brief, around 100-mg frozen tissues were minced and transferred to a 7-ml Dounce homogenizer (TIANDZ), to which was added 3 ml pre-chilled homogenization buffer, and kept on ice. After being allowed to stand for 5 − 10 min until full infiltration, tissues were vertically homogenized by the loose pestle (pestle A) until resistance was gone (around 10 − 15 times). The mixture was filtered using a 70-μm cell strainer into a pre-chilled 50-ml centrifuge tube. Two milliliter homogenization buffer was used to wash the homogenizer and cell strainer into the 50-ml centrifuge tube. After the homogenizer was cleaned using 2 ml nuclease-free H2O, the mixture in the 50-ml centrifuge tube was transferred into the cleaned homogenizer and further homogenized by 20 strokes of the tight pestle (pestle B). After this, the mixture was filtered through a 30-μm strainer into a 15-ml centrifuge tube and centrifuged at 500 g for 5 min at 4 °C to pellet nuclei. Pellets were resuspended using 1 ml blocking buffer into a 15-ml tube and centrifuged again at 500 g for 5 min at 4 °C. Nuclei pellets were then resuspended using 1.5 ml blocking buffer and mixed with a pipette. Nuclei suspensions were stained with 0.4% trypan blue to estimate whether nuclei were fully released under microscopic observation. Nuclei stained with DAPI were counted under microscope and samples with nuclei meeting the following criteria were used for library construction: concentration > 700/μl; volume ≥ 50 μl; membrane integrity rate ≥ 90%; agglomeration rate ≤ 5%; and clean solution with no visible impurities. The amygdala sample from the 6-year-old female rhesus macaque (SY1) was abandoned due to failure to meet the criteria, leaving a total of 39 brain samples suitable for further snRNA-seq data preparation (Fig. 1A, Additional file 2: Table S1).

Fig. 1
figure 1

Cellular census of rhesus macaque brain using snRNA-seq analysis. A Schematic overview of sample composition and anatomical positions of the ten macaque brain regions. F: female; M: male; AMY: amygdala; PU: putamen; HIP: hippocampus; TH: thalamus; DLPFC: dorsolateral prefrontal cortex; CG: cingulate gyrus; STG: superior temporal gyrus; SPL: superior parietal lobule; V4: visual cortex V4; CBC: cerebellar cortex. B Uniform Manifold Approximation and Projection (UMAP) plot of all 330,006 nuclei clustered into twenty major cell types. ExN: excitatory neuron; InN: inhibitory neuron; SPN: spiny projection neuron; CGC: cerebellar granule cell; TH-ExN: ExN mainly from TH; TH-InN: InN mainly from TH; AMY-ExN: ExN mainly from AMY; CA1–3: hippocampal excitatory CA1–3 cell; V4-ExN: ExN mainly from V4; Astro: astrocyte; Micro: microglia; OPC: oligodendrocyte progenitor cell; cOPC: differentiation committed OPC; Oligo: oligodendrocyte; Epen: ependymal cell; Endo: endothelial cell; Fib: perivascular fibroblast; SMC: smooth muscle cell; L5/6 NP: L5-6 near-projecting ExN; UL IT: upper-layer intratelencephalic-projecting ExN; DL CT: deep-layer corticothalamic-projecting ExN. C Dot plot shows expression of known markers or some new highly expressed genes of the twenty cell types. The dot size represents the percentage of nuclei expressing a gene by each cell type whereas dot color displays the scaled average expression level. CB-InN: InN mainly from CBC. D From left to right: dot plot displays significant annotation categories of cell type markers with color meaning − log-transformed adjusted P value and size meaning ratio of cell type markers in a term (first), while bar plots show number of cell types (second), percentage (perc.) of nuclei across individuals (third), regions (fourth) and age (fifth), and number of candidate cell type marker genes (sixth). ST: synaptic transmission; LCFA: long-chain fatty acid; dev.: development; GABA: gamma-aminobutyric acid; CNS: central nervous system; BBB: blood–brain barrier. E–M Validation of protein expression of cell-type-highly-expressed genes by immunofluorescence staining, including SV2B and CAMK2A for ExNs using slices from DLPFC (E); LMO4 and CAMK2A for ExNs using slices from SPL (F); PPP3CA and CAMK2A for ExNs using slices from SPL (G); PPP3CA and PPP1R1B for SPNs using slices from PU (H); CA12 and PPP1R1B for SPNs using slices from PU (I); PDE10A and PPP1R1B for SPNs using slices from PU (J); P2RY12 and RNASET2 for Micros using slices from white matter (K); OLIG2 and COL9A1 for OPCs using slices from STG (L); OLIG2 and VCAN for OPCs using slices from STG (M). Scale bars: 20 μm in E–J, K, and M; 10 μm in L

snRNA-seq library preparation

snRNA-seq libraries were prepared using DNBelab C Series Single-Cell Library Prep Set (MGI, 1000021082) as previously described [15, 27]. Briefly, following the manufacturer’s protocol, nuclei suspensions were successively used for droplet encapsulation, emulsion breakage, reverse transcription, and cDNA synthesis and amplification to generate barcoded libraries. The obtained libraries were then sequenced on the DNBSEQ platform at the China National GeneBank (Shenzhen, China) using the following sequencing strategy: Read 1 was 30 bp in length and included a 10-bp cell barcode 1, a 10-bp cell barcode 2, and a 10-bp unique molecular identifier (UMI). Read 2 was 90 bp in length for transcript sequence.

snRNA-seq data preprocessing

Raw DNBSEQ-sequencing reads were first filtered and demultiplexed using the parse tool in PISA (v0.2-17-g8cc39ef from https://github.com/MGI-tech-bioinformatics/DNBelab_C_Series_HT_scRNA-analysis-software) [28, 29]. Retained reads were aligned to the Macaca mulatta reference genome (Mmul_10 available at https://ftp.ensembl.org/pub/release-102/fasta/macaca_mulatta/dna/) [30] using STAR (v2.7.9a) [31], transferred to bam files with the sam2bam tool in PISA (v0.7 at https://github.com/shiquan/PISA) [29], and sorted by sambamba (v0.7.0) [32]. Gene annotation was performed for intron and exon reads using the anno tool in PISA (v0.7) with GTF file from Ensembl (v102 at https://ftp.ensembl.org/pub/release-102/gtf/macaca_mulatta/) [30]. UMIs were corrected using the corr tool in PISA (v0.2-17-g8cc39ef). Nucleus versus gene UMI count matrix was generated using the count tool in PISA (v0.7). Due to the issue of data format and environment configuration of the software, we used two PISA versions (i.e., v0.2-17-g8cc39ef and v0.7) with v0.7 being the primary one.

Nucleus clustering and cell-type annotation

First, each of the 39 samples was preprocessed using Seurat (v4.0.2) [33]. In brief, nuclei were filtered according to the following criteria: (1) genes were expressed by at least 100 nuclei; (2) a nucleus was required to express no fewer than 500 genes; (3) the number of UMIs in a nucleus was no fewer than 1000; (4) the percentage of mitochondrial genes expressed in a nucleus was not greater than 10%; and (5) the ratio of UMI counts to gene counts was no less than 1.2. Second, filtered datasets were dealt with using a series of Seurat functions, such as normalization with NormalizeData, data scaling with ScaleData, and dimension reduction with RunPCA. DoubletFinder [34] was employed to remove pseudo-doublets by assuming a doublet formation rate of 10%. Third, Seurat count objects were merged across all 39 brain samples and normalized using NormalizeData. FindVariableFeatures was performed to identify the 2000 most variable genes in terms of their expression. CellCycleScoring was applied to calculate S and G2-M scores based on the updated cell cycle markers [35]. ScaleData was used to scale the data and regress out the effects from sex, UMI counts, percentage of mitochondrial genes, and differences between the S and G2-M scores (i.e., “CC.differences”). We then performed RunPCA to reduce dimensions and RunHarmony to correct batch effects from individuals. Detection of Uniform manifold approximation and projection (UMAP) coordinates and clusters was performed using RunUMAP, FindNeighbors, and FindClusters sequentially. Cluster marker genes were identified using FindAllMarkers through a MAST test with UMI counts, CC.differences, and sex as covariates. We conducted FindAllMarkers using cells from all individuals as a whole analysis, or cells from each individual separately. A marker candidate was required to be expressed by at least 25% of the cluster cells (i.e., “min.pct”) and its Bonferroni-corrected P value and average fold change (FC) were < 0.05 and ≥ 1.5, respectively. We defined a gene as a cluster-specific marker if it was a marker candidate from both the whole analysis and separate analysis of at least two individuals. All cells were transferred to two published human and macaque brain snRNA-seq datasets [36, 37] using FindTransferAnchors and TransferData. For the published datasets, we kept cells as references which had anatomical origins same to our ten regions. One-to-one orthologous genes between human and rhesus macaque were identified using BioMart (Ensembl v102) [30]. Based on well-known marker genes and the two published NHP datasets, twenty major brain cell types were annotated. UMAP plots displaying different categories, such as cell types or regions, were performed using DimPlot. Functional annotation was performed for each cell type using their human coordinates of markers with the clusterProfiler (v4.0.5) R package [38] based on the “org.Hs.eg.db” dataset. We divided major cell types into nine groups: excitatory group including excitatory neurons (ExNs) mainly from the cortex, ExNs mainly from AMY (AMY-ExNs), hippocampal excitatory CA1–3 cells (CA1–3s) and ExNs mainly from V4 (V4-ExNs); inhibitory group including GABAergic inhibitory neurons (InNs); putamen-specific group including spiny projection neurons (SPNs); cerebellum-specific group including cerebellar granule cells (CGCs), InNs mainly from CBC (CB-InNs) and Bergmann glia (Bergmanns); thalamus-specific group including ExNs mainly from TH (TH-ExNs) and InNs mainly from TH (TH-InNs); astrocytic group including astrocytes (Astros); microglial group including microglia (Micros); oligodendrocyte-lineage group including oligodendrocyte progenitor cells (OPCs), differentiation committed OPCs (cOPCs) and oligodendrocytes (Oligos); membrane-related non-neuronal cell type group including ependymal cells (Epens), endothelial cells (Endos), perivascular fibroblasts (Fibs), and smooth muscle cells (SMCs). A similar pipeline was applied to identify the subtypes for each of these cell-type groups. Given few cells for some subtypes, we only performed FindAllMarkers using cells from all individuals. Two ways were applied to defined subtype markers: markers identified based on cells across subtypes from each cell type group were used to distinguish subtypes within the group, along with known markers of the cell type; markers identified based on cells across subtypes from all cell types were used to distinguish subtypes between cell-type groups.

Bulk RNA-seq data processing

The corresponding bulk RNA-seq raw data for the 39 snRNA-seq samples were obtained from our previous study [25]. The RNA-seq data were processed as previously described [39]. Briefly, we trimmed adaptors and low-quality reads using fastp (v0.21.0) [40] with the options “-n 15 -q 20 -u 40 -e 20 --length_required 40 -p -w 1.” Retained reads were aligned to the reference genome (Mmul_10) using HISAT2 (v2.2.1) [41] with parameters “-t --sensitive --no-discordant --no-mixed --–dta.” Expression of genes was estimated using StringTie2 (v2.1.4) with default parameters [42] based on the annotation file from Ensembl (v102). The sample versus gene matrix of fragments per kilobase of exon model per million mapped fragments (FPKM) values was utilized for the further correlation analysis.

Correlation of gene expression between snRNA-seq and bulk RNA-seq data

To assess consistency of gene expression between the snRNA-seq data and their corresponding bulk RNA-seq data, we performed pairwise correlation analysis between the two expression datasets using the cor.test command in R [43]. To construct the snRNA-seq expression matrix, UMI counts of each gene across all nuclei from a snRNA-seq sample were summed and transformed to FPKM values to represent the expression of the gene in that sample. Genes shared by the snRNA-seq and RNA-seq expression matrix were used. Both expression matrixes were log2 transformed by adding 1 prior to the correlation analysis. The bar plot as well as most plots of this study were plotted by the R package ggplot2 [44].

Percentage analysis

To compare cellular composition between distinct brain regions, age groups and individuals, we evaluated the percentages of nuclei per cell type or subtype across regions, age groups or individuals, respectively. The percentage of nuclei per region for a specific cell type was calculated by dividing the number of nuclei from a given cell type for a particular region by the total number of nuclei in that region. Similarly, the percentage of nuclei per individual for a specific subtype was calculated by dividing the number of nuclei from a given subtype for a specific individual by the total number of nuclei in that individual. When performing percentage analysis for subtypes across regions, we first normalized the number of nuclei per subtype from a region with the total number of nuclei from that region. We then calculated the percentage of nuclei per subtype for a specific region by dividing the normalized ratio per subtype from a given region by the total normalized ratios from all regions for that subtype. We performed proportion analysis for cell types or subtypes across age groups using a pipeline similar to that employed for regions. Cellular compositional difference of cell types was then calculated between the young and old groups using a linear regression model with sex as a covariate, setting the significant threshold to P value < 0.05.

Intercellular communication analysis

We independently conducted cell–cell communication analysis between cell types with regions and age being simultaneously considered or not, using the R package CellChat (v2.1.2) [45] by following the official workflow (https://github.com/sqjin/CellChat). For the ligand-receptor interaction database, we used “CellChatDB.human” that was manually curated in human and contained ~ 3300 validated molecular interactions involved in “paracrine/autocrine signaling,” “extracellular matrix (ECM)-receptor,” “cell–cell contact,” and “non-protein signaling.” Significant interactions were inferred through a permutation test on the probability value assigned to each interaction. We calculated the number of interactions between two combinations across regions and cell types, either considering age groups or not. We then estimated how numbers of intra- or inter-region cellular interactions correlated with regional spatial distance by linear regression using the R function lm [43]. The spatial distance value between two regions was measured following the widely used macaque brain atlas (http://www.brainmaps.org) [26]. The R function summary [43] was utilized to obtain the estimate regression coefficient (i.e., slope), R-squared (R2), and P values. The significant threshold of correlation was set to be P < 0.05. Slope values larger or less than zero represented positive or negative correlations, respectively.

Identification of region-specific genes (RSGs)

To identify RSGs, we first identified differentially expressed genes (DEGs) between nuclei from a region and all other nuclei from the remaining nine regions for each cell type using FindMarkers in Seurat through a MAST test, with thresholds of “min.pct” for the given cell type in the target region ≥ 0.25, Bonferroni-corrected P < 0.05, and an average FC ≥ 1.5. When conducting a MAST test using cells from all individuals as a whole analysis, UMI counts, CC.differences, and sex were treated as covariates. Genes defined as DEGs in multi-regions for the same cell type were discarded. The retained DEGs were considered as first-potential RSGs. In addition, we performed the above DEG analysis using cells from each individual separately, with UMI counts and CC.differences as covariates. We defined a first-potential RSG as a final-candidate RSG for a region per cell type if it was also a potential RSG for a same region and cell type in at least two individuals. Meanwhile, considering few or no cells from other regions for region-specific cell types (i.e., SPNs, CGCs, CB-InNs, Bergmanns, TH-ExNs, TH-InNs, AMY-ExNs, CA1–3s, and V4-ExNs), we regarded cell markers of these cell types as RSGs of the corresponding region and cell type. Functional annotation of RSGs was performed using the gprofiler2 R package [46].

Identification of aging-related DEGs

We identified DEGs between nuclei from the young and old macaques by cell type from each region through FindMarkers with a MAST test treating UMI, CC.differences, and sex as covariates. We first used cells from all individuals as a whole analysis. Then, we compared each of the two old individuals with each of the two young individuals to identify DEGs, namely SO1 vs. SY1, SO1 vs. SY2, SO2 vs. SY1, and SO2 vs. SY2. A gene was considered to be significantly upregulated in the old macaques if its “min.pct” of the given cell type in an old macaque region was ≥ 0.25, Bonferroni-corrected P value was < 0.05, and the average FC compared to the young macaques was ≥ 1.5. By contrast, a gene was considered to be significantly downregulated in the old macaques if its “min.pct” of the given cell type in a young macaque region was ≥ 0.25, Bonferroni-corrected P value was < 0.05, and the average FC compared to the young macaques was ≤ 1/1.5. An age-related DEG was required to exist both in the whole analysis and at least two of the four old versus young comparisons, with a same altering trend (namely simultaneously upregulated or downregulated) in a same region and cell type, but excluding those only from male-included comparisons (i.e., SO1 vs. SY2 and SO2 vs. SY2) to avoid sex-related effects, except for DEGs in AMY, as the young AMY sample was only from the young male individual. Then, we defined a DEG as a unique old-downregulated DEG if it was not upregulated in any cell type from any brain region of the old macaques. Conversely, we defined a DEG as a unique old-upregulated DEG if it was not downregulated in any cell type from any brain region of the old macaques. Functional enrichment analysis was performed using the gprofiler2 R package.

Transcriptional noise analysis

Transcriptional noise was evaluated following a previous pipeline [47]. Briefly, the numbers of nuclei from each cell type per region were down-sampled to be equal between young and old macaques and were no fewer than 10. Genes were divided into 10 equal bins based on their mean expression levels across the down-sampled nuclei with the top and bottom bins being excluded. For each bin, we only selected the 10% of genes with the lowest coefficient of variation (CV). We then calculated the euclidean distance between each cell and the corresponding mean value across nuclei from each cell type per region within each age group. This euclidean distance value was defined as transcriptional noise for each cell.

Pseudotime analysis

Pseudotime analysis was first performed for each cell type separately, using monocle3 [48, 49] by following the official tutorial (https://cole-trapnell-lab.github.io/monocle3/). Nuclei from different individuals were treated as distinct batches being corrected using the align_cds function. The root of a trajectory was programmatically or manually selected by using the node with enrichment of young cells. Genes varying over a trajectory or between two age groups were identified using the graph_test function with Moran’s I test by setting a significant threshold of adjusted P (Padj) < 0.05. If there were multiple trajectories split by subtypes for a cell type, we conducted a pseudotime analysis for each subtype, respectively. Expression dynamics of significantly varied genes along the trajectory were presented using the plot_genes_in_pseudotime function.

Heritability enrichment analysis

To ascertain correlations of complex human traits and diseases with specific signatures related to cell types, subtypes, regions, and aging, we independently performed linkage disequilibrium score regression (LDSC) (v1.0.1; https://github.com/bulik/ldsc) [50, 51] to predict the heritability enrichment of 30 diseases/traits across genes specifically or predominantly expressed in each cell type, subtype and region, or significantly differentially expressed between the young and old macaques as described previously [15, 52]. Cell-type- and subtype-specific marker genes, RSGs, and aging-related DEGs were transferred to hg19 genome coordinates based on orthologous genes from Ensembl. Files required for the basic baseline analysis were downloaded according to recommendations from LDSC authors (https://github.com/bulik/ldsc/wiki).

Immunofluorescence and immunohistochemistry staining experiments

Macaques for protein validation

Another two young and two aged rhesus macaques were used for immunohistochemistry and immunofluorescence validation. The two young macaques included SY3 (male, 5 years old) and SY4 (male, 11 years old), while the two old macaques included SO3 (female, 30 years old) and SO4 (male, 31 years old) (Additional file 2: Table S1). Macaque SY3 was used for validation of cell-type-specific genes; SY3 and SY4 were used for validation of most region-specific genes, except for ADCY5 where SO3 and SO4 were used due to the high expression of this gene in the old macaques; all four monkeys were used for validation of aging-related DEGs.

Brain acquisition, fixation and section

After euthanasia, each macaque brain was removed from the skull and immersed into 4% paraformaldehyde (PFA), 0.01 M PBS at 4 °C for 1 week. Then the 10 brain regions of interest from the right hemisphere were cut into 3 mm-thick sections coronally and embedded in paraffin. The paraffin-embedded sections were then sliced into 4-μm-thick slices on a Leica RM2016 Microtome (Leica, Shanghai, China) for immunohistochemistry and immunofluorescence staining.

Antibodies

The primary antibodies for immunohistochemistry and immunofluorescence staining included SV2B (SAB, Cat. #32139, 1:80), RNASET2 (SAB, Cat. #42738, 1:80), PDE10A (SAB, Cat. #44386, 1:100), PPP3CA (SAB, Cat. #32139, 1:80), VCAN (SAB, Cat. #43346, 1:80), COL9A1 (SAB, Cat. #46534, 1:80), ADCY5 (SAB, Cat. #34161, 1:200), HSP90AA1 (SAB, Cat. #32072, 1:100), P2RY12 (Cellsignal, Cat. #69766S, 1:400), CA12 (Cellsignal, Cat. #5865S, 1:100), ZIC4 (Abmart, Cat. #PH0841S, 1:200), OLIG2 (Abmart, Cat. #PA5829S, 1:200), TIAM1 (Solarbio, Cat. #K009405P, 1:200), ATP5MG (ZEN BIO, Cat. #825205, 1:40), ATP5MPL (ZEN BIO, Cat. #163326, 1:40), TCF7L2 (Abclonal, Cat. #A20770, 1:200), RGS16 (Abclonal, Cat. #A4078, 1:200), NTS (Abclonal, Cat. #A12326, 1:200), CFAP299 (Novus, Cat. #NBP1-86203, 1:200), CACNA1A (SAB, Cat. #37454, 1:200), HPCAL4 (SAB, Cat. #28970, 1:50), LMO4 (SAB, Cat. #39194, 1:80), VSNL1 (SAB, Cat. #47458, 1:50), PPP1R1B (SAB, Cat. #54693, 1:200) and CAMK2A (SAB, Cat. #46394, 1:200). The secondary antibodies for the immunohistochemistry or immunofluorescence staining contained anti-mouse/rabbit IgG detection system (ZSGB-BIO, Cat. #PV-9000), Cy2 AffiniPure donkey anti-rabbit IgG (Jackon ImmunoResearch, Cat. #715-225-150, 1:200), 488 AffiniPure goat anti-rabbit IgG (Jackon ImmunoResearch, Cat. #111-545-144, 1:200), Cy5 AffiniPure donkey anti-mouse IgG (Jackon ImmunoResearch, Cat. #715-175-150, 1:200) and Cy5 AffiniPure donkey anti-rabbit IgG (Jackon ImmunoResearch, Cat. #711-175-152, 1:200).

Immunohistochemistry staining

For the immunohistochemistry staining, paraffin sections were mounted on microscope slides and baked at 65 °C overnight. Then, the sections were deparaffinized and rehydrated with two 10-min changes of xylene and graded alcohols, respectively. The treated sections were later heated for 2 min with citrate buffer in a pressure cooker for antigen retrieval. Then the sections were treated with 70% formic acid for 7 min. Endogenous peroxidase in the sections was inactivated by incubation with 3% hydrogen peroxide in methanol for 10 min. The sections were then blocked with 2% bovine serum albumin (BSA) in 0.5% Triton X-100 for 30 min at 37 °C. The primary antibodies were diluted in blocking solution and incubated with the sections at 4 °C overnight. The sections were then incubated with secondary antibodies (PV-9000, 30 min) at 37 °C and developed using 3, 3’-diaminobenzidine (DAB, MXB, China) in chromogen solution, and counterstained with hematoxylin (Sigma, USA) for cell nucleus identification. In the negative-control experiments, only the primary antibodies were omitted. Finally, the sections were dehydrated and cleared with graded alcohols and xylene, respectively, and then covered with neutral gum and cover slips. The slides were analyzed using an Olympus CX41RF microscope (Olympus, Tokyo, Japan) and the images captured with the matching Olympus DP25 microscope digital camera. Immunostaining was performed in parallel in slides for each primary antibodies.

Immunofluorescence staining

The paraffin sections were mounted on microscope slides and baked at 65 °C overnight. Then, the sections were deparaffinized and rehydrated with two 10-min changes of xylene and graded alcohols, respectively. The treated sections were later heated for 2 min in citrate buffer in a pressure cooker for antigen retrieval. Then the sections were treated with 70% formic acid for 7 min. Endogenous peroxidase of the sections were inactivated by incubation with 3% hydrogen peroxide in methanol for 10 min. The sections were then blocked with 2% bovine serum albumin (BSA) in 0.5% Triton X-100 for 30 min at 37 °C. The primary antibodies were diluted in blocking solution and the sections were incubated at 4 °C overnight, followed by incubation with the corresponding fluorescent secondary antibodies at room temperature for 2 h. Nuclei were stained with 4’,6-diamidino-2-phenylindole (DAPI, Sigma, USA) for 15 min. In the negative-control experiments, the primary antibodies were omitted. The slides were mounted with Fluoromount-G (SouthernBiotech, USA) and imaged on a TissueFAXS Spectra System (TissueGnostics, Vienna, Austria). Immunostaining was performed in parallel in slides for each primary antibody.

Image analysis

Three coronal sections from each brain region were randomly selected for staining, and images of 3 visual fields were randomly collected for each section. All the images were captured with a × 20 objective under the same acquisition parameters for each immunostaining. The integrated optical density (IOD) of the immunohistochemistry staining and the proportion of immunopositive cells of the immunofluorescence staining in the 10 brain regions were analyzed quantitatively by ImageJ (NIH software). All images were converted into 8-bit gray scale and thresholds were adjusted and kept constant to highlight immunostained cells. Next, the immunoreactive particles were analyzed using the particle analysis tool. The IOD and ratio of the immunoreactive cells over each analyzed visual field were calculated. Values derived from a total of nine visual fields were averaged and presented as mean ± standard error of mean in bar plots.

Knockdown of VSNL1 and HPCAL4 in U251 glioma cells

To generate U251 glioma cell lines with knockdown of VSNL1 or HPCAL4, short hairpin RNA (shRNA) targeting specific genes were designed and cloned into the pLKO.1 lentiviral vector. The following target sequences were used for knockdown: sh-VSNL1-1: acagagtttaatgagcatgaa; sh-VSNL1-2: gcatgaactcaagcagtggta; sh-VSNL1-3: tgtccaagtgggaggctaaat; sh-HPCAL4-1: ccttgttcagaacactgagtt; sh-HPCAL4-2: tcagcagctctacatcaagtt; sh-HPCAL4-3: ccagattacattggaggagtt. Lentiviral particles were produced by co-transfecting the pLKO.1-shRNA plasmids along with packaging plasmids (psPAX2 and pMD2.G) into HEK-293 T cells using Lipofectamine 3000 (Thermo Fisher). The viral supernatants were collected 48 h post-transfection, filtered, and concentrated. U251 glioma cells were then transduced with the viral supernatants and selected using 2 μg/mL puromycin (InvivoGen, USA) for 3 days to ensure stable expression of the shRNAs. The efficiency of knockdown was assessed by quantitative PCR (qPCR), confirming the reduction of VSNL1 and HPCAL4 expression in transduced cells.

Senescence-associated beta-galactosidase (SA-β-gal) staining

SA-β-gal staining was used to identify senescent U251 glioma cells following the knockdown of VSNL1 or HPCAL4. This assay detects the enzymatic activity of SA-β-gal at a suboptimal pH of 6.0, a hallmark of senescent cells. In detail, U251 cells were transfected with shRNA targeting VSNL1 or HPCAL4 (or corresponding control shRNA), and after 48 h of culturing, they were subjected to SA-β-gal staining following the manufacturer’s protocol (C0602, Beyotime Biotech, Shanghai, China). The procedure was carried out as follows: (1) washing and fixation: cells were washed with hosphate-buffered saline (PBS) for 3 min, followed by fixation with the provided fixative solution at room temperature for 15 min; (2) staining: After fixation, cells were incubated with the SA-β-gal staining working solution overnight at 37 °C; (3) visualization: the stained cells were then visualized using the Cytation 5 cell-imaging multi-mode plate reader (BioTek, Winooski, VT, USA), which enabled the detection and quantification of SA-β-gal positive cells.

Quantitative real-time PCR (qPCR)

qPCR was performed to evaluate the knockdown efficiency of VSNL1 and HPCAL4 in glioma cells, as well as to measure the mRNA expression levels of senescence-related genes, including p16, p21, CCL2, CXCL1, CXCL3, and TNFα, in U251 glioma cells following the knockdown of VSNL1 or HPCAL4. Total RNA was extracted from U251 cells transfected with shRNA targeting VSNL1 or HPCAL4 (or control shRNA) using TRIzol reagent (15596018, Thermo Fisher Scientific, MA, USA) according to the manufacturer’s protocol. The RNA purity and concentration were assessed spectrophotometrically. For cDNA synthesis, an RT kit (K1622, Thermo Fisher Scientific, MA, USA) was used following the manufacturer’s instructions. qPCR was then performed to quantify the expression of target genes using specific primers for p16, p21, CCL2, CXCL1, CXCL3, TNFα, and β-actin (internal control), along with SYBR Green master mix (2 × Master qPCR Mix TSE201, Tsingke®, Beijing, China). The relative mRNA expression levels were determined using the 2−ΔΔCt method.

Statistical analyses

All statistical analyses were carried out based on data distribution. For immunofluorescence staining data of the eight RSGs, we performed the Anderson–Darling test on residuals for normality as the samples size was 60, using the R function ad.test (nortest) [53]. For immunohistochemical staining data of the eight aging-related genes and statistic data after knockdown of VSNL1 or HPCAL4 in U251 glioma cells, we separately performed the Shapiro–Wilk test for normality on raw values of each group or condition, using the R function shapiro.test [43]. Levene’s test was applied for checking homogeneity of variance, using the R function leveneTest (car) [54].

Consequently, data of RSGs such as RGS16, NTS, CA12, PDE10A, and TIAM1 were analyzed through one-way Analysis of Variance (ANOVA) followed by Dunnett multiple comparisons due to normal distribution and equal variances using the R functions of aov and glht (multcomp) [43, 55], while data of TCF7L2, ZIC4, and ADCY5 were analyzed through Welch’s ANOVA followed by Games-Howell test and adjusting P values with FDR due to normal distribution but unequal variances using the R functions of posthocTGH (userfriendlyscience) [56] and p.adjust [43]. Padj < 0.05 was set to be a significant threshold.

Data of aging-related genes such as ATP5MG, ATP5MPL (downregulated in old DLPFC), CFAP299, and VSNL1 were analyzed through unpaired two-tailed Student’s t test due to normal distribution and equal variances using the R function t.test [43], while data of CACNA1A and HSP90AA1 were analyzed through unpaired two-tailed Welch’s t-test due to normal distribution but unequal variances. Data of ATP5MPL (downregulated in old CG) and HPCAL4 were not normally distributed; thus, we analyzed through unpaired two-tailed Mann–Whitney U test using the R function wilcox.test [43]. The percentage of SA-β-gal-positive cells transduced with sh-VSNL1-1 plasmids, and relative qPCR expression levels of CXCL3 in the conditions of sh-VSNL1-1 and sh-VSNL1-3 were not normally distributed; thus, we compared their values with the condition of pLKO.1 using unpaired two-tailed Mann–Whitney U test, respectively. For data of other conditions from SA-β-gal staining and CXCL3 qPCR, the growth state of U251 cells (i.e., relative cell number), and relative qPCR expression levels of genes including p16, p21, CCL2, CXCL1, TNFα, VSNL1, and HPCAL4, unpaired two-tailed Student’s t test was used due to normal distribution and equal variances when compared with pLKO.1. P < 0.05 was set to be a significant threshold.

Results

Characterization of a cell atlas for NHP brains

To construct a single-cell transcriptome atlas for the NHP brain, we profiled gene expression in single nuclei of 39 samples from 10 regions of two young and two aged rhesus macaques using DNBelab C4 snRNA-seq (see “Methods” and Fig. 1A). After a series of filtration steps, a total of 330,006 high-quality nuclei were retained for further analyses, with a median of 2212 UMIs and 1398 genes per cell and of 6820 cells per sample (Additional file 1: Fig. S2). In terms of gene expression, the snRNA-seq data correlated well with their corresponding bulk RNA-seq data (Pearson correlation coefficient (PCC) = 0.71 ± 0.06; Additional file 1: Fig. S3 A), implying consistent transcriptomic profiles between snRNA-seq and bulk RNA-seq data. To correct for batch effects, the retained nuclei across all samples were integrated using Seurat with Harmony algorithm (Additional file 1: Fig. S3B–E). UMAP dimensional reduction and Louvain algorithm were applied to the integrated dataset to define distinct clusters, and DEGs were identified for each cluster based on the MAST test (Additional file 3: Table S2). Using canonical markers and cell-type transferring based on the published snRNA-seq datasets from human and macaque brains [36, 37], we identified twenty major cell types which were evenly split into neuronal and non-neuronal types (Fig. 1B–D, Additional file 1: Fig. S4). The neuronal cell types included excitatory neurons (ExNs; CAMK2A+ SLC17A7+) mainly from the cortex, GABAergic inhibitory neurons (InNs; GAD1+ GAD2+), ExNs mainly from V4 (V4-ExNs; POU6F2+ CUX2+), ExNs mainly from AMY (AMY-ExNs; DPYD+ PRR16+), ExNs mainly from TH (TH-ExNs; SLC17A6+ TCF7L2+), InNs mainly from TH (TH-InNs; GAD2+ KIT+), InNs mainly from CBC (CB-InNs; PVALB+ NEFH+), spiny projection neurons (SPNs; PPP1R1B+ PENK+) mainly from PU, cerebellar granule cells (CGCs; CBLN1+ NEUROD1+), and hippocampal excitatory CA1–3 cells (CA1–3s; CAMK2A+ LMO1+), while the non-neuronal types consisted of Bergmann glia (Bergmanns; FABP7+ NDRG2+), astrocytes (Astros; AGT+ AQP4+), microglia (Micros; CSF1R+ P2RY12+), oligodendrocyte progenitor cells (OPCs; OLIG2+ PCDH15+), differentiation committed OPCs (cOPCs; GPR17+ SOX4+), oligodendrocytes (Oligos; MOBP+ MOG+), ependymal cells (Epens; DNAH6+ FOXJ1+), endothelial cells (Endos; CLDN5+ FLT1+), perivascular fibroblasts (Fibs; SFRP2+ MYOC+), and smooth muscle cells (SMCs; ACTA2+ MYH11+). Functional annotation of the cell-type-specific genes revealed significant categories paralleling functions of the corresponding cell types (Fig. 1D). For example, ExN-specific genes were enriched in glutamatergic pathways, InN in GABA signaling pathways, SPN in dopamine response, Bergmann in cerebellar granular layer development, Astro in astrocyte differentiation, OPC, cOPC, and Oligo in gliogenesis and myelination, Micro in immune response, Epen in cilium assembly and cerebrospinal fluid circulation, and vascular cells in angiogenesis and wound healing. Besides, some cell types also showed intriguing functional enrichment, such as both TH-ExNs and TH-InNs, were remarkably associated with aerobic respiration, while TH-InNs were also specifically involved in the sterol and steroid metabolic processes. Astros were correlated to ossification and kidney development. Endos were related to maintenance of blood–brain barrier and aging.

In addition to the canonical cell-type markers, we also identified and validated many genes that were found to be predominantly expressed in a particular cell type but which, to our knowledge, have not yet been well reported in NHPs, including SV2B, PPP3CA, and LMO4 for neurons with an excitatory nature, PPP3CA also for SPNs, GRIK1 and GALNTL6 for InNs, CA12 and PDE10A for SPNs, HOMER3 and ZNF385C for CB-InNs, ANXA1 for Bergmanns, RGS16 for TH-ExNs, OTX2 and ASIC4 for TH-InNs, HPSE2 for V4-ExNs, ADGRV1 and ETNPPL for Astros, APBB1IP and RNASET2 for Micros, COL9A1 and VCAN for OPCs and cOPCs, LIMS2 and GNB4 for cOPCs, CD22 for cOPCs and Oligos, ANLN for Oligos, ARMC3 for Epens, and MFAP4 for Fibs (Fig. 1C, E–M, Additional file 3: Table S2). In particular, the expression of SV2B has been recently found to be restricted to the glutamatergic neurons of human temporal cortex [57]. Our results provide further evidence for SV2B being specifically expressed in primate ExNs. PPP3CA variants are known to be associated with severe epilepsy and dysmorphism [58]. Our finding of PPP3CA overexpression in the excitatory and spiny projection neurons of rhesus macaque might therefore provide a clue as to the pathogenic mechanism whereby these two cell types are involved in epileptic seizure. LMO4 has been shown to control the balance between excitatory and inhibitory neurons in the ventral spinal cord [59]. GRIK1 has recently been shown to be predominantly expressed in the GABAergic neurons in the rodent amygdala [60], while we confirmed its remarkable expression in macaque InNs. Local inactivation of GRIK1 expression in adult rodent amygdala reduces ongoing GABAergic transmission giving rise to a mild anxiety-like behavior [60]. ADGRV1 (alias VLGR1) variants have been reported to be related to epilepsy and audio-visual abnormalities [61]. High ADGRV1 expression in macaque astrocytes and astrocyte-like glial cells (i.e., Epens) would help to explain its pathological involvement. ETNPPL expression has been found to be specifically induced by diet in mouse brain astrocytes, while altered ETNPPL expression has been implicated in mood disorders [62, 63]. The predominant ETNPPL expression in macaque astrocytes may provide an insight into ETNPPL functions. Kelley et al. recently identified APBB1IP as a novel marker of human brain microglia even though its microglial functions remain unstudied [64]. Our own data show that APBB1IP expression is significantly associated with macaque microglia and support the view that APBB1IP may represent a useful microglial marker. Hamilton et al. have reported that microglial failure to digest apoptotic cells is a major contributor to RNASET2-deficient leukoencephalopathy [65]. RNASET2 overexpression in macaque Micros concurs with the view that RNASET2 pathogenesis is likely to operate through microglial dysfunction. Intriguingly, Pluvinage et al. found that CD22 is exclusively expressed by oligodendrocytes in human brain but by microglia in mouse brain [66]. Our own data on CD22 overexpression in macaque Oligos and cOPCs, a small population of oligodendroglial cells presenting an intermediate state between OPCs and newly formed Oligos, suggested that CD22 might be a specific marker of oligodendroglial cells in the primate brain. Consistent with our finding of high ANLN expression in macaque oligodendrocytes, Erwig et al. also stated that ANLN expression is at its highest in the myelinating oligodendrocytes of mouse brain and its depletion causes pathological myelin outfoldings by disturbing myelin septin assembly [67]. To our knowledge, functions of COL9A1 in the brain have not yet been reported, while its loss in mice disrupts myeloid cell functions [68]. Overexpression of COL9A1 in macaque OPCs and cOPCs should prompt future exploration of its functions. VCAN has been recently identified as a novel marker of a primary OPC subtype in a murine model of cerebral ischemia [69], but is also a marker of OPCs and an ExN subcluster in human brain [70]. Our results confirm that VCAN expression typically correlates with OPCs. Owing to the cost and precious nature of the samples involved, we opted to validate eight of these novel (or less well-known) cell-type-specific genes using immunostaining alongside canonical markers and obtained perfect concordance across macaque brain regions (Fig. 1E–M). In summary, our results expand the number of existing canonical markers for different brain cell types, including those specific to primates.

Cell–cell communications between brain cell types

Robust brain functions are critically reliant on the precise regulation of connections between cell types. Using CellChat, we analyzed intercellular communications across distinct cell types. In total, we observed 114 significant ligand-receptor pairs (LRPs) from 46 signaling pathways. The strongest crosstalk was evident among ExNs, InNs, CGCs, SPNs, and V4-ExNs (Fig. 2A). Cell types were found to be associated with the signaling pathways through five patterns, including two outgoing communication patterns of secreting cells and three incoming communication patterns of target cells (Fig. 2B). Four signaling groups were identified based on their functional similarity, as members within a group exhibited similar major senders and receivers and thus might perform similar and/or redundant roles (Fig. 2C–F). Group 1 was dominated by immune pathways (i.e., ADGRB, ANGPT, ANGPTL, BMP, CD39, CD45, COMPLEMENT, CX3C, ESAM, IL16, MHC-II, NTS, PECAM1, THY1, and VEGF) representing signaling from/to a single type of glial or vascular cells, especially Micros, Astros, and Endos. For example, the LRP ENTPD1-ADORA2B in the CD39 pathway specifically mediated signaling from Micros (Fig. 2F). ENTPD1 (alias CD39) is known as a microglial-specific marker [71], a characteristic that was also corroborated by our results (Fig. 2E). ENTPD1 encodes an ectonucleotidase that serves to regulate immunity and inflammation [72]. Both ESAM and PECAM1 exclusively regulate autocrine signaling between Endos (Fig. 2D), two endothelial cell markers that show selective expression in Endos from our data (Fig. 2E). ESAM encodes endothelial cell-selective adhesion molecule that regulates angiogenesis and endothelial permeability, with involvement in various vascular diseases including atherosclerosis and various hypervascular tumors [73, 74]. PECAM1 plays a vital role in the maintenance of endothelial junction integrity and the vascular barrier [75]. Group 2 included twelve pathways (i.e., CLDN, EPHA, JAM, LAMININ, MPZ, NEGR, Netrin, PDGF, PTN, PTPR, RELN, SLIT), which were mainly involved in nervous system development and largely represented signaling from neuronal or non-neuronal cells as primary outgoing resources. For instance, NEGR1 was particularly evident in terms of the connections among neurons (Fig. 2D, F). NEGR1 encodes neuronal growth regulator 1, a cell adhesion molecule that is involved in cortical layering and which is essential for the balance between excitatory/inhibitory neurons [76]. NEGR1 deficiency dramatically alters brain morphology and disrupts neurite arborization in mice [76]. Moreover, the LRP PTN-PTPRZ1 in the PTN pathway primarily output signals from non-neuronal cells, consistent with the very low level of PTN expression in neurons (Fig. 2D–F). PTN encodes pleiotrophin, a growth factor that modulates microglia-mediated neuroinflammation and is involved in many brain disorders including neurodegenerative diseases and drug addiction [77]. Group 3 contained nine pathways (i.e., ADGRL, CADM, CNTN, NCAM, NRG, NRXN, PCDH, PSAP, and PTPRM) with predominant signals from neuronal cells and cOPCs, which are mainly associated with cell–cell adhesion, synapse assembly, and behavior. In particular, the LRP NRXN3-NLGN1 in the NRXN pathway primarily sends signals from neuronal cells to most cell types (Fig. 2F). NRXN3 was highly expressed in neuronal cells and one of its missense variant has been shown to enhance mouse empathy fear [78]. NLGN1 encodes a postsynaptic cell adhesion protein that influences synapse formation and stabilization and has been linked to cognitive disorders, such as schizophrenia, autism, and AD [79, 80]. Group 4 contained mainly integrin-mediated signaling pathways (i.e., ADGRE, CDH, COLLAGEN, FN1, GAP, IGF, PROS, NOTCH, TENASCIN) targeting to Astros and SMCs, which have been implicated in cell junction assembly and wound healing. For example, both LRPs FN1-SDC4 in the FN1 pathway and TNR-SDC4 in the TENASCIN pathway used SDC4 as a receptor and send signals to Astros (Fig. 2F). SDC4 and FN1 were identified as a marker specific to Astros and vascular cells in our data, respectively (Fig. 2E), both of which play important roles in cell adhesion, cell migration, and wound healing [81]. Together, the intercellular communication analysis improved our understanding of how brain cells connect with each other to maintain normal brain functions.

Fig. 2
figure 2

Cell–cell communications across the twenty major cell types. A Interaction weights/strength between any pair of cell types, displayed by edge width. B Alluvial plots highlight the outgoing signaling patterns of secreting cells (left) and incoming signaling patterns of target cells (right). C Signaling pathways were clustered into four groups based on their functional similarity. Dot size was proportional to the overall communication probability (Prob.). D Outgoing communication patterns of secreting cells (left) and incoming communication patterns of target cells (right). Dots were colored according to cell types, while dot size was proportion to the contribution score (Contri.) of each cell type to each signaling pathway to show association between cell type and their enriched signaling pathways. Signaling pathways were colored according to the four functionally similar groups. E Dot plot illustrates expression of the representative ligand-receptor genes for each of the four functionally similar groups. The dot size represents the percentage of nuclei expressing a gene by each cell type whereas dot color displays the scaled average expression level. F Interaction weights/strength between any pair of cell types for the representative ligand-receptor pairs, displayed by edge width

Transcriptional landscape and regional specificity of cell subtypes

Vast cellular heterogeneity within the neuronal population has been previously reported in both primate and mouse brain [82, 83]. Revealing the transcriptomic diversity of brain cell types is a prerequisite for the exploration of brain functions that are dependent upon cellular and regional specificity. Thus, we performed further cluster analysis for the major cell types so as to define their subtypes using a similar method. In total, we detected 88 subclusters, including 3 for SPNs, 2 for CB-InNs (i.e., Purkinje neurons and basket/stellate cells), 6 for TH-ExNs, 2 for TH-InNs, 2 for AMY-ExNs, 4 for CA1–3s, 5 for V4-ExNs, 12 for cortical ExNs, 15 for InNs, 5 for Astros, 9 for Micros, 2 for OPCs, 6 for Oligos, 2 for Epens, 10 for Endos, and 3 for Fibs (Fig. 1B, Additional file 1: Figs. S5–13). These subpopulations of cells were characterized by high expression of their own specific signatures (Additional file 3: Table S2). Thus, for example, we detected excitatory subclasses mainly from cortical neurons based on known layer markers (i.e., CUX2 for L2- 3, RORB for L3-5, RXFP1, TLE4 and FOXP2 for L5/6, and HTR2C for L5-6 near-projecting (NP) ExNs) and annotations from the published human snRNA-seq dataset (Additional file 1: Fig. S4) [37]. Interestingly, most V4 L2-5 ExNs were distinct from excitatory layer neurons of cortex tissues like DLPFC, CG, STG, and SPL, while some hippocampal CA ExNs were close to cortical ExNs. InNs were roughly split into two major branches based on their developmental origins in the medial and caudal ganglionic eminence using LHX6 and ADARB2 as a typical marker, respectively. LHX6 exhibits high chromatin at the promoter site of SST and PVALB, while ADARB2 had accessibility at the promoter chromatin of VIP and LAMP5, which are four well-known inhibitory subtype markers [27] (Additional file 1: Fig. S9). Besides, specific signatures were identified between the two cerebellar inhibitory subtypes (Additional file 1: Fig. S6), including their well-known markers (i.e., PCP2, CALB1, and CA8) and unwell-characterized genes (e.g., PPP1R17, GRIK1, and COL4 A2) for Purkinje neurons, and unwell-characterized genes (e.g., KIT, BTBD11, and TFAP2B) for basket/stellate cells which have been often defined using its well-known markers shared with Purkinje neurons. However, we found specific over-expression of PCP2, CALB1, and CA8 in macaque Purkinje neurons rather than basket/stellate cells, while both subtypes highly expressed CB-InN markers, such as GAD1, PVALB, NEFH, HOMER3, and ZNF385C. Vascular and ependymal cells are two major membrane-related non-neurogenic cell types. Based on known markers and the public human brain vasculature dataset [84], we confirmed the three well-defined vascular populations including Endos, Fibs, and the mural cell SMCs (Additional file 1: Fig. S13). Endos were grouped into 10 subtypes, including 7 capillary clusters (capEndos by MFSD2A+) as well as one arteriole cluster (aEndos by VEGFC+ ARL15+ DKK2+ FBLN5+), one venule cluster (vEndos by ADGRG6+ ACKR1+ DTX4+), and one new Endos (KHDRBS2+ CAMK2A+ CABP1+). Fibroblasts highly expressed MFAP4 in rhesus macaques and were classed into three clusters, including Fib1 (CEMIP+ KCNK2+ SLC6A20+), Fib2 (SFRP2+ MYOC+ GABRA6+ GRM4+), and Fib3 (SFRP2+ MYOC+ TP63+ ADAMTSL1+). Ependymal cells were clustered into Epen1 (FAM183A+ LHB+ C1H1orf87+) and Epen2 (HTR2C+ PRLR+ MAP3K15+).

Intriguingly, we observed some region-specific subtypes (RSSs) that consistently overexpressed a series of genes across distinct cell types, especially with PU, TH, and CBC (Fig. 3A, Additional file 1: Figs. S5–13D). For example, SPN1–3 and capEndo4 were primarily from PU and simultaneously expressed CA12, PDE10A, PENK, PPP1R1B, and ADCY5; TH-ExN1–6, Astro4, Micro5, and Oligo5 were mainly from TH and highly expressed NTS, RGS16, SLC17A6, and TCF7L2; Purkinje neurons, basket/stellate cells, CGCs, Bergmanns, Micro6, Oligo6, and Fib2 were predominantly from CBC and widely expressed CHN2, GRID2, PATJ, TIAM1, and ZNF521. We also found that some RSSs exhibited highly distinctive transcriptomic signatures. Thus, TH-specific InN1 specifically expressed ASIC4 and OTX2. We further validated the predominant protein expression of these genes in corresponding brain regions, such as PDE10A in PU, NTS in TH, and TIAM1 in CBC (Fig. 3B–E, Additional file 1: Figs. S14 A and S15). PDE10A encodes a cyclic adenosine monophosphate (cAMP)/cyclic guanosine monophosphate (cGMP) phosphodiesterase and is associated with dopaminergic neurotransmission in the striatum (comprising PU and caudate nucleus), while its dysfunction causes various neurological and psychiatric disorders such as AD, PD, and schizophrenia [85]. NTS encodes a neuropeptide, neurotensin, that has been found to be highly expressed in mouse TH and is essential for sleep regulation through a thalamo-amygdala circuit in mice [86]. OTX2 encodes a homeodomain transcription factor that plays a crucial role in the identity and fate of thalamic glutamatergic precursors by inhibiting GABAergic differentiation [87]. TIAM1 encodes T-cell lymphoma invasion and metastasis-inducing protein 1, a Rac-specific activator critical for the migration of cerebellar granule cell precursors along a gradient of brain-derived neurotrophic factor (BDNF) [88]. The functions of these genes suggest that the region-specific cell subtypes may contribute to the functional specificity of brain regions and brain disorders.

Fig. 3
figure 3

Characterization of region-specific subtypes (RSSs). A Heat map (top) indicates percentage of nuclei counts across the ten brain regions for each cell subtype, whereas the dot plot (bottom) shows expression of RSS marker genes. The dot size represents the percentage of nuclei expressing a gene by subtype whereas dot color displays the average expression level. B–E Immunofluorescence validation of protein expression of RSS-overexpressed genes, including NTS for TH (B and C) and TIAM1 for CBC (D and E) across the ten regions. Bar plots in C and E represent the mean percentage of immunopositive cells for each staining across three visual fields with standard error of mean bars. For both genes, one-way analysis of variance (ANOVA) was applied followed by Dunnett multiple comparisons due to normal distribution and equal variances. ****: Padj < 0.0001; ***: Padj < 0.001; **: Padj < 0.01; *: Padj < 0.05. Scale bars: 20 μm

Delineation of region-specific expression patterns shared among different cell types

The aforementioned signatures shared across RSSs inform region-specific expression and promise to provide new insights into the molecular mechanisms underlying the regionalization of brain functions. To further explore brain region-specific expression patterns, we first evaluated the cellular composition of different regions and found that brain cortex tissues (e.g., STG, SPL, V4, CG, DLPFC, and CBC) harbored a little higher proportion of neurons than subcortical nuclei (e.g., TH, HIP, PU, and AMY) whereas the latter contained more glial cells (Fig. 4A, B). Moreover, subcortical neurons were distinct from cortical neurons and exhibited their own signatures, particularly remarkable expression of known layer markers in cortical excitatory neurons but not in subcortical excitatory neurons. These results corroborate the observation that neuron addition paralleled primate cortical expansion, especially in the generation of upper-layer neurons [89, 90], whereas glia are the chief orchestrators of neurodegenerative diseases related to subcortical nuclei, such as glial involvement in PD pathogenesis within basal ganglia, an important component of which is putamen [91]. In addition, most cerebellar cells were also very different from other brain tissue cells and thus named separately, such as the excitatory neurons—CGCs, inhibitory neurons—Purkinje and basket/stellate cells, and astrocytes—Bergmanns, which indicates the unique roles of CBC.

Fig. 4
figure 4

Cellular and molecular patterns across the ten regions. A Percentage of cell types in each brain region. B Percentage of neuron, glia, and vascular in cortical and subcortical tissues. C Correlations between regional expression Euclidean distance and spatial distance by linear regression. P and R2 values were obtained using the “lm” and “summary” functions in R. D Scatter plot of slope and R2 values from “lm” to show correlations between intra- or inter-region cellular communication and spatial distance. Dot size represents P values. E Number of region-specific genes (RSGs) in each region by totals (top) or shared by the 1–20 cell types (bottom). The bottom dot plot indicates the number of RSGs by dot size; its y-axis represents the number of cell types in which an RSG was simultaneously identified. F Dot plot shows expression of top AMY-, PU-, TH-, V4-, and CBC-specific genes across regions. The dot size represents the percentage of nuclei expressing a gene by each region whereas dot color displays the scaled average expression level. G–J Immunofluorescence validation of RSG protein expression, including CA12 for PU (G and H) and TCF7L2 for TH (I and J) across the ten regions. Bar plots in H and J represent the mean percentage of immunopositive cells for each staining across three visual fields with standard error of mean bars. One-way ANOVA was applied for CA12 followed by Dunnett multiple comparisons due to normal distribution and equal variances, while Welch’s ANOVA was applied for TCF7L2 followed by Games-Howell test and adjusting P values using FDR due to normal distribution but unequal variances. ****: Padj < 0.0001; ***: Padj < 0.001; **: Padj < 0.01; *: Padj < 0.05. Scale bars: 20 μm

To explore how transcriptomic profiles changed with spatial distance variation among brain regions, we first estimated expression Euclidean distance between two regions across genes. We detected significantly positive correlations (slope > 0, P = 0.02066, R2 = 0.03128, linear model) between regional expression and spatial distance (Fig. 4C, Additional file 1: Fig. S16 A), which suggests that more distant regions exhibited greater transcriptomic differences. Moreover, we assessed intra- or inter-region cellular communications using CellChat and calculated relationships between the numbers of LRPs and spatial distance through a linear model. We found that connections between CBC-specific (i.e., CGCs, CB-InNs and Bergmanns) and other cell types were extensively significantly increased with increasing spatial distance, but connections within CBC-specific cell types tended to be the opposite (Fig. 4D, Additional file 1: Fig. S17). Meanwhile, connections with glia and other neurons mainly exhibit negative relationships with spatial distance. These imply that cerebellar cells had strong interactions with cells from those distant regions, while other cellular interactions between regions tended to decrease with increasing spatial distance. The differences in gene expression and cell–cell communications between regions may make a significant contribution to the functional regionalization of brain.

Thus, to further characterize the regional transcriptomic specificity, we identified region-specific genes (RSGs) that were significantly highly expressed in a cell type of one region as compared to the other nine regions. In total, we obtained 1532 unique RSGs, with most of them deriving from a particular cell type and also being defined as an RSG in a same region from previously published NHP datasets [36, 92, 93] (Fig. 4E, Additional file 1: Fig. S16B, Additional file 4: Table S3). TH, CBC, AMY, and PU had the largest number of RSGs. Enrichment analysis of RSGs revealed fundamental biological pathways that were common among several regions, including nervous system development, cell communication, cognition, and behavior (Additional file 1: Fig. S16 C). In addition, we also annotated unique categories relevant to regional functions. For example, RSGs in AMY regulate axon development and response to temperature stimulus, while AMY is a brain region primarily involved into emotional response [94]; RSGs in PU are related to visual learning and behavior, while PU has been found to be implicated in visual information process [95]; RSGs in HIP are associated with long-term memory; RSGs in TH are involved in cellular respiration and the adenosine triphosphate (ATP) metabolic process, whereas TH is a high energy-consuming gateway that relays information flow to the cortex [96]; RSGs in V4 play important roles in response to retinoic acid and short-term memory. Interestingly, we noted that most RSS-shared overexpressed genes were also cell-type-shared RSGs in a same region, including NTS, RGS16, and TCF7L2 in TH, ADCY5, CA12, and PDE10A in PU, and TIAM1 in CBC, whose proteins were also highly expressed in the corresponding brain regions (Figs. 3 and 4F–J, Additional file 1: Fig. S14). RGS16 encodes a member of the G protein signaling (RGS) family and is predominantly expressed in the rodent thalamus and suprachiasmatic nucleus. RGS16 knockdown causes abnormal circadian locomotor activity rhythms and food anticipation [97]. TCF7L2, encoding a member of the LEF1/TCF transcription factor family, modulates postmitotic differentiation and terminal selection postnatally in TH and is a risk gene for mental disorders like schizophrenia and autism [98]. ADCY5 encodes adenylate cyclase 5, a striatum-specific enzyme that converts ATP to cAMP. Mutations in ADCY5 are associated with dyskinesias, while PU and striatum are engaged in the modulation of movement and learning [99]. CA12 encodes membrane-associated carbonic anhydrase XII which catalyzes the reversible hydration of carbon dioxide; its expression is affected by hypoxia and estrogen receptors [100]. Despite its still unclear functions in PU, CA12 is regarded as a cancer marker and a potential biomarker for chronic back pain, while increased left PU volume is associated with the chronic back pain caused by ankylosing spondylitis [100, 101]. We also validated highly protein expression of ZIC4 in macaque CBC, a CBC-specific gene being defined in human [92, 93]. ZIC4 encodes a zinc-finger (ZF) transcription factor closely related to the protein product of the CBC-specific gene ZIC1. Both ZIC1 and ZIC4 are required for cerebellar development and their heterozygous deletion causes reduced cerebellar size due to decreased proliferation of postnatal granule cell progenitors and is associated with Dandy-Walker malformation (DWM), the most common congenital cerebellar defect [102]. Overall, these results are a reflection of the specific molecular and cellular signatures that exist among different cell types in the same brain region; these may play important roles in maintaining the regional transcriptional network thereby contributing to the functional integrity and regionalization of the brain.

Cellular alterations in the aging brain of NHPs

Aging is the progressive decline in physiological and functional capacity, involving especially memory loss, decreased comprehension and a slowing of reflexes [14, 103, 104]. Previous studies have used single-cell/nucleus RNA sequencing to show that aging is driven by various complex molecular and cellular signatures from a few brain regions [47, 105, 106]. Here, using a wider range of brain regions, we set out to survey the transcriptional and cellular changes that accompany aging. We first determined the correlations of inter-region expression and spatial distance as described above for the young and old macaques, respectively. We found significant positive relationships between regional expression and spatial distance in the young group but not in the old individuals (Fig. 5A), implying that differences in gene expression between regions largely decreased in the old macaques. We then compared the transcriptional noise between young and old individuals and detected broadly increased instability of gene expression with aging (Additional file 1: Fig. S18), consistent with the reduced transcriptional stability apparent in aged tissues that has been reported previously [47, 107].

Fig. 5
figure 5

Alterations in transcriptomic profiles of macaque brain with aging at single-cell levels. A Correlations between regional expression Euclidean distance and spatial distance by linear regression for the young and old macaques, respectively. P and R2 values were obtained using the “lm” and “summary” functions in R. B Circos plots show the number of DEGs downregulated (left panel) and upregulated (right panel) in the old macaques compared to the young macaques across different cell types and brain regions. Each connecting curve denotes the number of intersected down- or upregulated DEGs between two cell-type-region groups. The connecting curves suggest similar altering transcriptomic trends of aging-associated DEGs across regions and cell types in an age group. Venn diagram (middle panel) indicates the number of unique or shared old-downregulated or old-upregulated DEGs. C Comparison of protein expression for selected DEGs in corresponding young and old brain regions by immunohistochemistry staining. The old-downregulated DEGs included ATP5MG in DLPFC; ATP5MPL in DLPFC and CG; HSP90AA1, VSNL1, and HPCAL4 in CG; and CFAP299 in CBC. CACNA1A were old-upregulated in CBC. Bar plots (bottom panels) show the mean integrated optical density (IOD) for each immunohistochemistry staining across nine visual fields with standard error of mean bars. Scale bars: 50 μm. D The growth state of U251 glioma cells transduced with pLKO.1, sh-VSNL1-2, and sh-HPCAL4-2 plasmids, respectively. Red arrows pointed cells that increased in size. The right bar plot shows the mean relative cell number for each treatment with standard error of mean bars. E The senescence-associated beta-galactosidase (SA-β-gal) positive state of U251 glioma cells transduced with pLKO.1, sh-VSNL1-3, and sh-HPCAL4-3 plasmids, respectively. The right bar plot shows the mean percentage of SA-β-gal-positive cells for each treatment with standard error of mean bars. F Relative mRNA expression levels of senescence marker genes (i.e., p16, p21, CCL2, CXCL1, CXCL3, and TNFα) based on qPCR in U251 glioma cells after each treatment of VSNL1 or HPCAL4 knockdown. In C–F, unpaired Welch’s t-test was applied for CACNA1A and HSP90AA1 due to normal distribution but unequal variances and unpaired Mann–Whitney U test was applied for ATP5MPL (CG), HPCAL4, SA-β-gal staining in sh-VSNL1-1, and CXCL3 qPCR expression in sh-VSNL1-1 and sh-VSNL1-3 due to abnormal distribution, while other comparisons were all analyzed using unpaired Student’s t test due to normal distribution and equal variances. ****: P < 0.0001; ***: < 0.001; **: P < 0.01; *: P < 0.05; ns: P > 0.05

We independently performed intra- or inter-region cell–cell communication analysis for the young and old groups (Additional file 1: Fig. S19). The total interaction strength of the inferred cell–cell communication networks from the old samples was higher than those from the young samples (Additional file 1: Fig. S19 A). Compared to the young group, interaction strength of most cell types from different regions was enhanced in the old group, while some other types were attenuated (Additional file 1: Fig. S19B). For example, Oligos from CBC and CG, and CA1–3s from HIP showed increased outgoing and incoming signaling in the old samples, while Endos from CBC, SMCs from STG, and Micros from V4 exhibited the opposite. Signaling pathways highly enriched in the old macaques were a little more than those in the young individuals (Additional file 1: Fig. S19 C), with the top ones involved in signal transduction and pathogenesis of cancer and neurodegenerative diseases such as L1CAM, OPIOID, WNT, CSF, KIT, CCK, THBS, TRAIL, and ncWNT. The positive relationship of spatial distance with connections between CBC-specific and most other cell types were detected in both young and old group but enhanced in the latter (Additional file 1: Fig. S19D). Most cross-regional connections to TH-InNs, InNs, or Astros were significantly attenuated with increasing spatial distance in the young macaques, but the pattern eliminated in the old macaques. These results suggest that cross-regional cell–cell communication largely changed in aging macaques, especially cellular connections in CBC to other regions. So far, cerebellum has exhibited tight links with aging and neurodegenerative disorders [52, 108, 109]. Our findings may therefore provide an important new perspective on understanding the molecular underpinnings of the aging process and age-associated pathologies.

To estimate aging difference in cell-type composition, we then performed differential proportion analysis between the young and old macaques using a linear regression model with sex as a covariate. We found a large increase of CGCs and AMY-ExNs as well as decrease of InNs, cOPCs, and Oligos in the old group (Additional file 1: Fig. S19E). Meanwhile, several subtypes were also obviously either reduced or enriched in the old group (Additional file 1: Figs. S5–13E). For example, although no remarkable difference was found in the whole composition of Micros and OPCs between the young and old groups, their subtypes Micro8 and OPC1 significantly decreased in the old macaques. These results suggest that inhibitory neurons and glia were largely reduced with aging, while excitatory neurons in CBC and AMY were the opposite. The cellular compositional changes might contribute to the fundamental aging mechanisms. For instance, decline in inhibitory neurons and glia might correlate with age-related reduced functional plasticity and immune dysfunction, respectively, as dendritic remodeling of inhibitory neurons and immune response of immunocompetent glial cells diminishes with aging [110, 111].

To explore the potential molecular mechanisms underlying brain aging, we examined the gene expression profile alterations associated with aging across brain regions and cell types. In total, we identified 1386 and 1276 genes that were respectively significantly down- and upregulated in the old macaques compared to the young macaques (Fig. 5B, Additional file 5: Table S4). These aging-related DEGs exhibited similar trends of expression changes across regions and cell types in an age group, with 138 DEGs simultaneously down- and upregulated in the old macaques at different cell types or regions. Many of these DEGs were included as human aging-related genes by the Aging Atlas database [112] or work of Jia et al. based on large-scale Genotype-Tissue Expression project (GTEx) samples [113] (Additional file 1: Fig. S20 A). Age-related transcriptomic changes were largely found in ExNs, InNs, Astros, and Oligos (Fig. 5B), with marginally more old-downregulated DEGs in InNs, Micros, OPCs, CGCs, and Bergmanns than their old-upregulated DEGs. Compared to the other eight regions, DLPFC and CG showed the most remarkable changes in gene expression with aging, in line with other evidence for the important roles played by the two regions in aging, such as the involvement of DLPFC dysfunction in the cognitive deficits associated with age-related disorders [114] and abnormal neuronal and glial activity of CG in the pathogenesis of AD depression [115]. Further, we separately conducted functional enrichment analysis using DEGs from each cell type and region. Genes displaying decreased expression in the old macaques across cell types and regions were extensively enriched in energy metabolism, apoptosis, translation, immunity, and gliogenesis, whereas DEGs upregulated in the aging macaques were more predominantly implicated in signal transduction, synapse assembly, dendrite development, neurogenesis, nervous system process, cognition, behavior, sensory perception of pain, and rhythm (Additional file 1: Fig. S20B). Nevertheless, opposite exceptions also widely existed. For instance, upregulated DEGs in SPNs from PU and Oligos from TH were implicated in gliogenesis, while downregulated DEGs from many cell types and regions also showed relatively weak enrichment of cell communication and signal transduction. Genes involved in immunity were mostly downregulated in glial and vascular cells (especially Astros) of the old macaques, with some exceptions upregulated in the old, such as Micros in SPL and V4, Oligos in TH, and Endos in CG. These results suggest heterogeneous aging phenotypes in distinct cell types across different brain regions. Moreover, we observed significant expression changes in several genes with aging that were then confirmed by immunohistochemical staining showing that these changes were reflected at the protein level; these included ATP5MG and ATP5MPL significantly downregulated in the DLPFC and CG of old macaques, HSP90AA1, VSNL1, and HPCAL4 downregulated in the old CG, CFAP299 downregulated in the old CBC as well as CACNA1A upregulated in the old CBC (Fig. 5C, Additional file 1: Figs. S20 C and S21). Both ATP5MG (alias ATP5L) and ATP5MPL (alias ATP5MJ) are potential regulators of ATP synthesis and function in energy metabolism [116, 117]. Interestingly, expression of ATP5MG is also decreased in AD, a well-known aging-associated disease [116]. HSP90AA1 encodes a heat shock protein that participate in autophagy or immunity, and which has potential therapeutic or prognostic value in the context of disease [118, 119]. Both VSNL1 (alias VILIP-1) and HPCAL4 (alias VILIP-2) are members of the visinin-like subfamily of neuronal calcium sensor proteins [120, 121]. VSNL1 encodes visinin-like protein 1, a peripheral AD biomarker. VSNL1 upregulates cAMP levels and differentiation in peripheral cells and is associated with schizophrenia [122]. Very little is known about the function of HPCAL4, which may affect neurotransmitter release through modulating the presynaptic P/Q-type Ca2 + channels [120]. We identified HPCAL4 as an ExN marker and VSNL1 as a marker of ExNs, TH-ExNs, and TH-InNs using all cells from the ten brain regions (Additional file 1: Fig. S22 A, Additional file 3: Table S2), while both genes were identified as ExN markers based on only cells from CG (Additional file 1: Fig. S22B). We also validated their protein expression in ExNs from CG and TH (only for VSNL1) using immunofluorescence staining (Additional file 1: Fig. S22 C). These results imply potential important links between functions of VSNL1 and HPCAL4 and neurons, especially ExNs. Few functional data are available for CFAP299 (alias C4orf22), which is mainly known for its possible involvement in spermatogenesis through regulating cell cycle and apoptosis [123]. CACNA1A encodes the pore-forming subunit of neuronal P/Q-type Ca2 + channels and its mutations may cause cognitive impairment, autism, and epileptic encephalopathy [124].

To examine effects of aging-related DEGs on aging, we conducted loss-of-function experiments in U251 glioma cells for VSNL1 and HPCAL4, respectively. We found that knockdown of VSNL1 or HPCAL4 both significantly induced U251 cells senescence (Fig. 5D–F, Additional file 1: Figs. S23 and S24): cell proliferation remarkably slowed down; cell volume obviously increased; percentage of senescence-associated beta-galactosidase (SA-β-gal) positive cells largely enhanced; expression of senescence marker genes (i.e., p16, p21, CCL2, CXCL1, CXCL3, and TNFα) significantly increased. These results suggest that VSNL1 and HPCAL4 downregulation in the old macaques might contribute to aging.

Moreover, to further investigate cellular transcriptional alterations with aging, we performed pseudotime analysis to infer the aging process for a cell type or subtype using monocle3 (Additional file 1: Fig. S25). Many cell types or subtypes exhibited trajectory transition from young to aged nuclei, which suggest a gradual progression toward aging in these cell types/subtypes, such as Astros, SPNs in PU, and cerebellum-specific cell types/subtypes like CGCs, Bergmanns, and Purkinje neurons. However, not all cells from the young individuals retained a youthful state, while a subset of cells from the old individuals located in an earlier pseudotime. Nevertheless, many aging-related DEGs showed a similar expression trend across pseudotime. For instance, the old-downregulated DEGs like CFAP299 and ATP5MPL exhibited a decrease expression along pseudotime, while the old-upregulated DEGs like CACNA1A had an increase expression along pseudotime. These results indicate that a progressive aging trajectory potentially existed in most cells of many cell types.

Taken together, our results provide an unparalleled close-up view of the landscape of transcriptomic and cellular alterations that accompany aging, especially female aging due to our sample composition as we discussed detailedly in the following “Discussion” section, in an NHP brain.

Implications of cellular, regional, and aging-relevant signatures in human diseases and traits

To assess the potential involvement of gene expression signatures related to specific cell types, subtypes, regions, and aging in human diseases and complex traits, we used LDSC to determine the heritability enrichment of 30 diseases/traits across genes that were specifically or predominantly expressed in each cell type, subtype, and region, or which exhibited differential expression after aging, as previously described [15, 52] (Fig. 6, Additional file 1: Fig. S26). As was to be expected, neurological traits were substantially enriched in neuronal populations, including worry, sleep duration, and schizophrenia primarily in excitatory and spiny projection neurons, intelligence, and years of education in excitatory and inhibitory neurons, as well as insomnia and neuroticism mainly in InNs. We also observed enrichment of schizophrenia in astrocytes, consistent with the observation that astrocyte dysfunction contributes critically to the pathogenesis of schizophrenia [125]. In addition, we also detected enriched heritability of immune-related diseases, such as rheumatoid arthritis and type 1 diabetes, in several microglial subtypes. In particular, enrichment of attention-deficit/hyperactivity disorder (ADHD) was mainly found in V4-ExNs and Micros. ADHD is a common neurodevelopmental disorder in children and adults, characterized by inattention, impulsivity, and hyperactivity. Visual cortex abnormalities have been detected in adults with ADHD, while alterations in microglial activation have been implicated in the ADHD pathophysiology [126, 127]. Our results further confirmed these findings and revealed potential involvements of V4-ExNs and Micros in ADHD pathogenesis. Among the three lipid-related traits, i.e., levels of high-density lipoprotein (HDL), low-density lipoprotein (LDL), and triglycerides, we found enrichment of HDL in Astros but not for the other two traits, in line with what is known about the regulation of astrocytes in HDL generation [128]. Moreover, we also observed enrichment of neurological traits across aging-related genes in regions including DLPFC, CG, and V4, implying that genes which exhibit significant expression changes with aging may play important roles in the development of human diseases and complex traits, especially in the context of mental illness. Triglycerides and coronary artery disease (CAD) showed highest heritability enrichment in DLPFC glia across aging-related genes, suggesting potential roles of DLPFC glia in abnormal triglyceride levels and increased CAD risk with aging.

Fig. 6
figure 6

Heritability enrichment of 30 human complex traits and diseases across specific signatures of cell types and subtypes. Heat maps show the association of selected human traits and diseases with cell types (left of the top panel) and subtypes (right of the top panel) using − log-transformed P values from LDSC. Circle size indicates significance with empty space denoting P ≥ 0.05. The rectangles colored by deep sky blue highlight the selected patterns of enrichment

Taken together, neuronal cells were more likely to be involved in neurologically relevant traits, with distinct subtypes contributing to different traits. Microglial markers were engaged in immune-related traits. Aging-related DEGs across cell types from several regions were substantially involved in neurological traits, suggesting that the expression of many genes implicated in brain function and brain disorders may change during the aging process. By contrast, no obvious enrichment patterns of the traits were found in the brain region-specific genes, while aging-related genes did not appear to be involved in traits that are generally regarded as being systemic, such as height and LDL. In summary, the 30 diseases/traits analyzed here exhibited tight connections to cell-type markers and aging-related DEGs rather than region-specific genes in the macaque brain, especially for neurological or immune-related traits.

Discussion

The molecular and cellular diversity and specializations in human and mouse brains have been extensively studied by scRNA-seq or snRNA-seq, significantly advancing our understanding of the organization and dynamics of brain networks underlying cognition and behavior [19, 37, 129]. Employing snRNA-seq, we systematically profiled the cellular and molecular architecture of the rhesus macaque brain across multiple distinct regions. Our findings offer a comprehensive delineation of the cellular and molecular organization of NHP brain regions, providing new insights into how a broad functional repertoire emerges from static brain anatomy—an enduring challenge in brain research [129].

Through differential gene expression analysis and immunofluorescence, we expanded the list of markers for major brain cell types and identified several subtypes with distinctive marker genes. These cell populations were not isolated but engaged in intricate communication via diverse LRPs. Additionally, we uncovered cellular and molecular mechanisms underlying brain regionalization. For instance, differences in gene expression across brain regions increased with spatial distance between regions. Notably, cerebellar cells exhibited enhanced interactions with cells in distant brain regions, whereas interactions among most other cell types tended to decrease with distance. A comparison of cellular composition revealed a higher proportion of neurons in the cortical regions compared to subcortical nuclei, where glial cells were more abundant. Moreover, we observed several cell subtypes that exhibited regional specificity and identified a number of genes that were highly expressed in different cell types within the same region. The region-specific subtypes commonly overexpressed many region-specific genes. Thus, we envisage that region-specific subtypes and genes may play crucial roles in conferring functional and structural specificity to brain networks across regions.

Furthermore, we observed significant transcriptional and cellular alterations in the brain with aging. Compared to the young macaques, aged macaques exhibited reduced transcriptional stability and decreased gene expression differences across regions. Cross-regional cell–cell communication also showed marked changes with aging, particularly in interactions involving cerebellar cells. Certain cell types displayed aging-related compositional alterations, including an increase in CGCs and AMY-ExNs, and a decrease in InNs, cOPCs, and Oligos. Thousands of aging-related DEGs were identified across regions and cell types, with extensively consistent trends in their expression changes. Functional enrichment analysis revealed that upregulated DEGs in aging macaques were associated with cognition and behavior, while downregulated DEGs were primarily involved in energy metabolism and apoptosis. Pseudotime analysis showed that most cell types in aged macaques followed a progressive aging trajectory, with many age-related DEGs exhibiting consistent expression patterns over pseudotime. These results highlight extensive molecular and cellular changes as critical drivers of aging, characterized by weakened regional specialization and connections, disrupted functional properties, and impaired cognition and behavior. Such changes might be implicated in age-related diseases like AD, PD, and schizophrenia [13, 14, 116, 125]. Finally, heritability enrichment analyses of human complex traits and diseases across cell populations, brain regions, and age groups revealed a strong association of neurological traits with neuronal populations and aging-related DEGs, particularly in regions such as DLPFC, CG, and V4. Our dataset serves as a valuable resource for investigating normal brain function and the dysfunctions associated with aging and mental disorders. However, our samples were mainly from female individuals, and to avoid sex-related effects, we identified age-related DEGs by incorporating multiple comparisons with those only from male-to-female compares being discarded, which suggested that the aging-related alterations we reported here primarily reflected female aging. More samples from both male and female individuals in the further would make up for this limitation. Even more comprehensive surveys will now be required to determine how the alterations in transcriptomic signatures that occur with age are involved in the process of aging and aging-associated pathogenesis. For instance, knockdown of VSNL1 or HPCAL4 in U251 cells induced remarkable senescence, suggesting that the downregulation of these genes in aged macaques may contribute to aging. However, further experimental and bioinformatic analyses are required to clarify their precise roles in the aging process and associated pathologies.

Conclusions

It has always been fascinating and challenging to investigate how brain functions with precision and has evolved rapidly as the most complex and vital organ of humans. Yet, due to ethical constraints and limited access to human brain tissues, NHPs that present genetic, physiological, and anatomical similarities to humans make ideal models for studying brain function and evolution. By characterizing the cellular transcriptome profiles of the rhesus macaque across multiple regions from young and aging individuals using snRNA-seq, this work unveils the intricate molecular and cellular organization of the primate brain and its dynamic alterations with aging, especially in females. The data generated in this study would therefore provide a valuable resource for the in-depth analysis of the mechanisms underlying the precise modulation of brain function, the evolution of human high cognitive abilities, and the progression of brain aging and related diseases.

Data availability

Processed data have been deposited into the Science Data Bank database at https://doiorg.publicaciones.saludcastillayleon.es/https://doiorg.publicaciones.saludcastillayleon.es/10.57760/sciencedb.06819 [130]. Raw sequence data have been deposited in the Genome Sequence Archive under accession number CRA009154 at https://bigd.big.ac.cn/gsa/browse/CRA009154 [131, 132]. The code for the primary analyses in this manuscript is available at https://github.com/WanYMEN/Regional-and-aging-specific-cellular-architecture-of-non-human-primate-brains [133]. The bulk RNA-seq raw data from Li et al. [25] are available at https://ngdc.cncb.ac.cn/gsa/browse/CRA000336. The snRNA-seq dataset of macaque brain from Chiou et al. [36] are available at https://cellxgene.cziscience.com/collections/8c4bcf0d-b4df-45c7-888c-74fb0013e9e7. The snRNA-seq dataset of human brain from Siletti et al. [37] are available at https://cellxgene.cziscience.com/collections/283d65eb-dd53-496d-adb7-7570c7caa443. Region-specific genes from Kang et al. [92] and the BrainBase database [93] are available at https://www.nature.com/articles/nature10523 and https://ngdc.cncb.ac.cn/brainbase/index, respectively. Human aging-related genes from the Aging Atlas database [112] and Jia et al. [113] are available at https://ngdc.cncb.ac.cn/aging/age_related_genes and https://www.nature.com/articles/s41420-018-0093-y, respectively.

Abbreviations

AD:

Alzheimer’s disease

ADHD:

Attention-deficit/hyperactivity disorder

aEndo:

Arteriole endothelial cell

AMY:

Amygdala

AMY-ExN:

ExN mainly from AMY

Astro:

Astrocyte

CA1–3:

Hippocampal excitatory CA1–3 cell

CAD:

Coronary artery disease

capEndo:

Capillary endothelial cell

CB-InN:

InN mainly from CBC

CBC:

Cerebellar cortex

CG:

Cingulate gyrus

CGC:

Cerebellar granule cell

cOPC:

Differentiation committed OPC

CV:

Coefficient of variation

DEG:

Differentially expressed gene

DL CT:

Deep-layer corticothalamic-projecting ExN

DLPFC:

Dorsolateral prefrontal cortex

Endo:

Endothelial cell

Epen:

Ependymal cell

ExN:

Excitatory neuron

FC:

Fold change

Fib:

Perivascular fibroblast

FPKM:

Fragments per kilobase of exon model per million mapped fragments

HDL:

High-density lipoprotein

HIP:

Hippocampus

IBD:

Inflammatory bowel disease

InN:

Inhibitory neuron

IOD:

Integrated optical density

L5/6 NP:

L5-6 near-projecting ExN

LDL:

Low-density lipoprotein

LDSC:

Linkage disequilibrium score regression

LRP:

Ligand-receptor pair

Micro:

Microglia

NHP:

Non-human primate

Oligo:

Oligodendrocyte

OPC:

Oligodendrocyte progenitor cell

P adj :

Adjusted P value

PCC:

Pearson correlation coefficient

PD:

Parkinson’s disease

PU:

Putamen

qPCR:

Quantitative PCR

RSG:

Region-specific gene

RSS:

Region-specific subtype

SA-β-gal:

Senescence-associated beta-galactosidase

scRNA-seq:

Single-cell RNA sequencing

shRNA:

Short hairpin RNA

SMC:

Smooth muscle cell

snRNA-seq:

Single-nucleus RNA sequencing

SPL:

Superior parietal lobule

SPN:

Spiny projection neuron

STG:

Superior temporal gyrus

TH:

Thalamus

TH-ExN:

ExN mainly from TH

TH-InN:

InN mainly from TH

UL IT:

Upper-layer intratelencephalic-projecting ExN

UMAP:

Uniform manifold approximation and projection

UMI:

Unique molecular identifier

V4:

Visual cortex V4

V4-ExN:

ExN mainly from V4

vEndo:

Venule endothelial cell

References

  1. Rogers J, Gibbs RA. Comparative primate genomics: emerging patterns of genome content and dynamics. Nat Rev Genet. 2014;15(5):347–59.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Ma Q, Ma W, Song T-Z, Wu Z, Liu Z, Hu Z, et al. Single-nucleus transcriptomic profiling of multiple organs in a rhesus macaque model of SARS-CoV-2 infection. Zool Res. 2022;43(6):1041–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Fan C, Wu Y, Rui X, Yang Y, Ling C, Liu S, et al. Animal models for COVID-19: advances, gaps and perspectives. Signal Transduct Target Ther. 2022;7(1):220.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Krienen FM, Goldman M, Zhang Q, C. H. del Rosario R, Florio M, Machold R, et al. Innovations present in the primate interneuron repertoire. Nature. 2020;586(7828):262–269.

  5. Khrameeva E, Kurochkin I, Han D, Guijarro P, Kanton S, Santel M, et al. Single-cell-resolution transcriptome map of human, chimpanzee, bonobo, and macaque brains. Genome Res. 2020;30(5):776–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Passingham R. How good is the macaque monkey model of the human brain? Curr Opin Neurobiol. 2009;19(1):6–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Dolan RJ. Emotion, cognition, and behavior. Science. 2002;298(5596):1191–4.

    Article  CAS  PubMed  Google Scholar 

  8. Sjöstedt E, Zhong W, Fagerberg L, Karlsson M, Mitsios N, Adori C, et al. An atlas of the protein-coding genes in the human, pig, and mouse brain. Science. 2020;367(6482):eaay5947.

  9. Poo M-m. Transcriptome, connectome and neuromodulation of the primate brain. Cell. 2022;185(15):2636–9.

    Article  CAS  PubMed  Google Scholar 

  10. Velmeshev D, Schirmer L, Jung D, Haeussler M, Perez Y, Mayer S, et al. Single-cell genomics identifies cell type–specific molecular changes in autism. Science. 2019;364(6441):685–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Notaras M, Lodhi A, Dündar F, Collier P, Sayles NM, Tilgner H, et al. Schizophrenia is defined by cell-specific neuropathology and multiple neurodevelopmental mechanisms in patient-derived cerebral organoids. Mol Psychiatry. 2022;27(3):1416–34.

    Article  CAS  PubMed  Google Scholar 

  12. Mathys H, Davila-Velderrain J, Peng Z, Gao F, Mohammadi S, Young JZ, et al. Single-cell transcriptomic analysis of Alzheimer’s disease. Nature. 2019;570(7761):332–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Hindle JV. Ageing, neurodegeneration and Parkinson’s disease. Age Ageing. 2010;39(2):156–61.

    Article  PubMed  Google Scholar 

  14. Xia X, Jiang Q, McDermott J, Han J-DJ. Aging and Alzheimer’s disease: Comparison and associations from molecular to system level. Aging Cell. 2018;17(5):e12802.

  15. Han L, Wei X, Liu C, Volpe G, Zhuang Z, Zou X, et al. Cell transcriptomic atlas of the non-human primate Macaca fascicularis. Nature. 2022;604(7907):723–31.

    Article  CAS  PubMed  Google Scholar 

  16. Consortium* TS, Jones RC, Karkanias J, Krasnow MA, Pisco AO, Quake SR, et al. The Tabula Sapiens: A multiple-organ, single-cell transcriptomic atlas of humans. Science. 2022;376(6594):eabl4896.

  17. Schaum N, Karkanias J, Neff NF, May AP, Quake SR, Wyss-Coray T, et al. Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature. 2018;562(7727):367–72.

    Article  PubMed Central  Google Scholar 

  18. Yuan Y, Sun D-M, Qin T, Mao S-Y, Zhu W-Y, Yin Y-Y, et al. Single-cell transcriptomic landscape of the sheep rumen provides insights into physiological programming development and adaptation of digestive strategies. Zool Res. 2022;43(4):634–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Saunders A, Macosko EZ, Wysoker A, Goldman M, Krienen FM, de Rivera H, et al. Molecular diversity and specializations among the cells of the adult mouse brain. Cell. 2018;174(4):1015-1030.e16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Wei J-R, Hao Z-Z, Xu C, Huang M, Tang L, Xu N, et al. Identification of visual cortex cell types and species differences using single-cell RNA sequencing. Nat Commun. 2022;13(1):6902.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Hao Z-Z, Wei J-R, Xiao D, Liu R, Xu N, Tang L, et al. Single-cell transcriptomics of adult macaque hippocampus reveals neural precursor cell populations. Nat Neurosci. 2022;25(6):805–17.

    Article  CAS  PubMed  Google Scholar 

  22. Zhou X, Lu Y, Zhao F, Dong J, Ma W, Zhong S, et al. Deciphering the spatial-temporal transcriptional landscape of human hypothalamus development. Cell Stem Cell. 2022;29(2):328-343.e5.

    Article  CAS  PubMed  Google Scholar 

  23. Ziffra RS, Kim CN, Ross JM, Wilfert A, Turner TN, Haeussler M, et al. Single-cell epigenomics reveals mechanisms of human cortical development. Nature. 2021;598(7879):205–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Habib N, Li Y, Heidenreich M, Swiech L, Avraham-Davidi I, Trombetta JJ, et al. Div-Seq: Single-nucleus RNA-Seq reveals dynamics of rare adult newborn neurons. Science. 2016;353(6302):925–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Li M-L, Wu S-H, Zhang J-J, Tian H-Y, Shao Y, Wang Z-B, et al. 547 transcriptomes from 44 brain areas reveal features of the aging brain in non-human primates. Genome Biol. 2019;20(1):258.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Mikula S, Trotts I, Stone JM, Jones EG. Internet-enabled high-resolution brain mapping and virtual microscopy. Neuroimage. 2007;35(1):9–15.

    Article  PubMed  Google Scholar 

  27. Lei Y, Cheng M, Li Z, Zhuang Z, Wu L, sun Y, et al. Spatially resolved gene regulatory and disease-related vulnerability map of the adult macaque cortex. Nat Commun. 2022;13(1):6747.

  28. Liu C, Wu T, Fan F, Liu Y, Wu L, Junkin M, et al. A portable and cost-effective microfluidic system for massively parallel single-cell transcriptome profiling. bioRxiv. 2019:818450.

  29. Shi Q, Liu S, Kristiansen K, Liu L. The FASTQ+ format and PISA. Bioinformatics. 2022;38(19):4639–42.

    Article  CAS  PubMed  Google Scholar 

  30. Yates AD, Achuthan P, Akanni W, Allen J, Allen J, Alvarez-Jarreta J, et al. Ensembl 2020. Nucleic Acids Res. 2019;48(D1):D682–8.

    PubMed Central  Google Scholar 

  31. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.

    Article  CAS  PubMed  Google Scholar 

  32. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: fast processing of NGS alignment formats. Bioinformatics. 2015;31(12):2032–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573-3587.e29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. McGinnis CS, Murrow LM, Gartner ZJ. DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 2019;8(4):329-337.e4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Tirosh I, Izar B, Prakadan SM, Wadsworth MH, Treacy D, Trombetta JJ, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352(6282):189–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Chiou KL, Huang X, Bohlen MO, Tremblay S, DeCasien AR, O’Day DR, et al. A single-cell multi-omic atlas spanning the adult rhesus macaque brain. Science Advances. 9(41):eadh1914.

  37. Siletti K, Hodge R, Mossi Albiach A, Lee KW, Ding S-L, Hu L, et al. Transcriptomic diversity of cell types across the adult human brain. Science. 382(6667):eadd7046.

  38. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3).

  39. Wang Y-M, Ye L-Q, Wang M-S, Zhang J-J, Khederzadeh S, Irwin DM, et al. Unveiling the functional and evolutionary landscape of RNA editing in chicken using genomics and transcriptomics. Zool Res. 2022;43(6):1011.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Kovaka S, Zimin AV, Pertea GM, Razaghi R, Salzberg SL, Pertea M. Transcriptome assembly from long-read RNA-seq alignments with StringTie2. Genome Biol. 2019;20(1):278.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Team RC. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2022. https://www.R-project.org/.

  44. Wickham H. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York; 2016. https://ggplot2.tidyverse.org.

  45. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan C-H, et al. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12(1):1088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Raudvere U, Kolberg L, Kuzmin I, Arak T, Adler P, Peterson H, et al. g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 2019;47(W1):W191–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Zhang H, Li J, Ren J, Sun S, Ma S, Zhang W, et al. Single-nucleus transcriptomic landscape of primate hippocampal aging. Protein Cell. 2021;12(9):695–716.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32(4):381–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Cao J, Spielmann M, Qiu X, Huang X, Ibrahim DM, Hill AJ, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature. 2019;566(7745):496–502.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Bulik-Sullivan BK, Loh PR, Finucane HK, Ripke S, Yang J, Patterson N, et al. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat Genet. 2015;47(3):291–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Finucane HK, Bulik-Sullivan B, Gusev A, Trynka G, Reshef Y, Loh P-R, et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat Genet. 2015;47(11):1228–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Finucane HK, Reshef YA, Anttila V, Slowikowski K, Gusev A, Byrnes A, et al. Heritability enrichment of specifically expressed genes identifies disease-relevant tissues and cell types. Nat Genet. 2018;50(4):621–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Gross J, Ligges U. nortest: Tests for Normality. R package version 1.0–4; 2015. https://CRAN.R-project.org/package=nortest.

  54. Fox J, Weisberg S. An {R} Companion to Applied Regression, Third Edition. Thousand Oaks CA: Sage; 2019. https://r-forge.r-project.org/projects/car.

  55. Hothorn T, Bretz F, Westfall P. Simultaneous Inference in General Parametric Models. Biom J. 2008;50(3):346–63.

    Article  PubMed  Google Scholar 

  56. Peters G-J. Userfriendlyscience: Quantitative analysis made accessible. R package version 0.7.2; 2018. https://doiorg.publicaciones.saludcastillayleon.es/10.17605/osf.io/txequ.

  57. Pazarlar BA, Aripaka SS, Petukhov V, Pinborg L, Khodosevich K, Mikkelsen JD. Expression profile of synaptic vesicle glycoprotein 2A, B, and C paralogues in temporal neocortex tissue from patients with temporal lobe epilepsy (TLE). Mol Brain. 2022;15(1):45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Rydzanicz M, Wachowska M, Cook EC, Lisowski P, Kuźniewska B, Szymańska K, et al. Novel calcineurin A (PPP3CA) variant associated with epilepsy, constitutive enzyme activation and downregulation of protein expression. Eur J Hum Genet. 2019;27(1):61–9.

    Article  CAS  PubMed  Google Scholar 

  59. Joshi K, Lee S, Lee B, Lee JW, Lee S-K. LMO4 Controls the Balance between Excitatory and Inhibitory Spinal V2 Interneurons. Neuron. 2009;61(6):839–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Englund J, Haikonen J, Shteinikov V, Amarilla SP, Atanasova T, Shintyapina A, et al. Downregulation of kainate receptors regulating GABAergic transmission in amygdala after early life stress is associated with anxiety-like behavior in rodents. Transl Psychiatry. 2021;11(1):538.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Zhou P, Meng H, Liang X, Lei X, Zhang J, Bian W, et al. ADGRV1 variants in febrile seizures/epilepsy with antecedent febrile seizures and their associations with audio-visual abnormalities. Front Mol Neurosci. 2022;15: 864074.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. White CJ, Ellis JM, Wolfgang MJ. The role of ethanolamine phosphate phospholyase in regulation of astrocyte lipid homeostasis. J Biol Chem. 2021;297(1): 100830.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Shao L, Vawter MP. Shared gene expression alterations in schizophrenia and bipolar disorder. Biol Psychiatry. 2008;64(2):89–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Kelley KW, Nakao-Inoue H, Molofsky AV, Oldham MC. Variation among intact tissue samples reveals the core transcriptional features of human CNS cell classes. Nat Neurosci. 2018;21(9):1171–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Hamilton N, Rutherford HA, Petts JJ, Isles HM, Weber T, Henneke M, et al. The failure of microglia to digest developmental apoptotic cells contributes to the pathology of RNASET2-deficient leukoencephalopathy. Glia. 2020;68(7):1531–45.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Pluvinage JV, Sun J, Claes C, Flynn RA, Haney MS, Iram T, et al. The CD22-IGF2R interaction is a therapeutic target for microglial lysosome dysfunction in Niemann-Pick type C. Sci Transl Med. 2021;13(622):eabg2919.

  67. Erwig MS, Patzig J, Steyer AM, Dibaj P, Heilmann M, Heilmann I, et al. Anillin facilitates septin assembly to prevent pathological outfoldings of central nervous system myelin. eLife. 2019;8:e43888.

  68. Probst K, Stermann J, Bomhard I, Etich J, Pitzler L, Niehoff A, et al. Depletion of collagen IX alpha1 impairs myeloid cell function. Stem Cells. 2018;36(11):1752–63.

    Article  CAS  PubMed  Google Scholar 

  69. Shen S, Wang Y, Xia R, Zeng F, Cao J, Liu X, et al. Differentiation independent neuroprotective role of Vcan+ oligodendrocyte precursor cells in poststroke cognitive impairment recovery. bioRxiv. 2020:2020.11.07.356311.

  70. Tran MN, Maynard KR, Spangler A, Huuki LA, Montgomery KD, Sadashivaiah V, et al. Single-nucleus transcriptome analysis reveals cell-type-specific molecular signatures across reward circuitry in the human brain. Neuron. 2021;109(19):3088-3103.e5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Butovsky O, Weiner HL. Microglial signatures and their role in health and disease. Nat Rev Neurosci. 2018;19(10):622–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Antonioli L, Pacher P, Vizi ES, Haskó G. CD39 and CD73 in immunity and inflammation. Trends Mol Med. 2013;19(6):355–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Inoue M, Ishida T, Yasuda T, Toh R, Hara T, Cangara HM, et al. Endothelial cell-selective adhesion molecule modulates atherosclerosis through plaque angiogenesis and monocyte–endothelial interaction. Microvasc Res. 2010;80(2):179–87.

    Article  CAS  PubMed  Google Scholar 

  74. Cangara HM, Ishida T, Hara T, Sun L, Toh R, Rikitake Y, et al. Role of endothelial cell-selective adhesion molecule in hematogeneous metastasis. Microvasc Res. 2010;80(1):133–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Privratsky JR, Newman PJ. PECAM-1: regulator of endothelial junctional integrity. Cell Tissue Res. 2014;355(3):607–19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Singh K, Jayaram M, Kaare M, Leidmaa E, Jagomäe T, Heinla I, et al. Neural cell adhesion molecule Negr1 deficiency in mouse results in structural brain endophenotypes and behavioral deviations related to psychiatric disorders. Sci Rep. 2019;9(1):5457.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Fernández-Calle R, Vicente-Rodríguez M, Gramage E, Pita J, Pérez-García C, Ferrer-Alcón M, et al. Pleiotrophin regulates microglia-mediated neuroinflammation. J Neuroinflammation. 2017;14(1):46.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Keum S, Kim A, Shin JJ, Kim J-H, Park J, Shin H-S. A Missense Variant at the Nrxn3 Locus Enhances Empathy Fear in the Mouse. Neuron. 2018;98(3):588-601.e5.

    Article  CAS  PubMed  Google Scholar 

  79. Dufort-Gervais J, Provost C, Charbonneau L, Norris CM, Calon F, Mongrain V, et al. Neuroligin-1 is altered in the hippocampus of Alzheimer’s disease patients and mouse models, and modulates the toxicity of amyloid-beta oligomers. Sci Rep. 2020;10(1):6956.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Qin L, Liu Z, Guo S, Han Y, Wang X, Ren W, et al. Astrocytic Neuroligin-3 influences gene expression and social behavior, but is dispensable for synapse number. Mol Psychiatry. 2024.

  81. Maldonado H, Savage BD, Barker HR, May U, Vähätupa M, Badiani RK, et al. Systemically administered wound-homing peptide accelerates wound healing by modulating syndecan-4 function. Nat Commun. 2023;14(1):8069.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Ma S, Skarica M, Li Q, Xu C, Risgaard RD, Tebbenkamp ATN, et al. Molecular and cellular evolution of the primate dorsolateral prefrontal cortex. Science. 2022;377(6614):eabo7257.

  83. Yao Z, van Velthoven CTJ, Nguyen TN, Goldy J, Sedeno-Cortes AE, Baftizadeh F, et al. A taxonomy of transcriptomic cell types across the isocortex and hippocampal formation. Cell. 2021;184(12):3222-3241.e26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Garcia FJ, Sun N, Lee H, Godlewski B, Mathys H, Galani K, et al. Single-cell dissection of the human brain vasculature. Nature. 2022;603(7903):893–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Amin HS, Parikh PK, Ghate MD. Medicinal chemistry strategies for the development of phosphodiesterase 10A (PDE10A) inhibitors - An update of recent progress. Eur J Med Chem. 2021;214: 113155.

    Article  CAS  PubMed  Google Scholar 

  86. Ma C, Zhong P, Liu D, Barger ZK, Zhou L, Chang W-C, et al. Sleep regulation by neurotensinergic neurons in a thalamo-amygdala circuit. Neuron. 2019;103(2):323-334.e7.

    Article  PubMed  Google Scholar 

  87. Puelles E, Acampora D, Gogoi R, Tuorto F, Papalia A, Guillemot F, et al. Otx2 controls identity and fate of glutamatergic progenitors of the thalamus by repressing GABAergic differentiation. J Neurosci. 2006;26(22):5955.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Zhou P, Porcionatto M, Pilapil M, Chen Y, Choi Y, Tolias KF, et al. Polarized signaling endosomes coordinate BDNF-induced chemotaxis of cerebellar precursors. Neuron. 2007;55(1):53–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  89. Elston GN, Benavides-Piccione R, Elston A, Zietsch B, Defelipe J, Manger P, et al. Specializations of the granular prefrontal cortex of primates: Implications for cognitive processing. Anat Rec A Discov Mol Cell Evol Biol. 2006;288A(1):26–35.

    Article  Google Scholar 

  90. Molyneaux BJ, Arlotta P, Menezes JRL, Macklis JD. Neuronal subtype specification in the cerebral cortex. Nat Rev Neurosci. 2007;8(6):427–37.

    Article  CAS  PubMed  Google Scholar 

  91. Bhaduri B, Alladi PA. Glial cells as key orchestrators of neural degeneration in basal ganglia disorders. In: Patro I, Seth P, Patro N, Tandon PN, editors. The Biology of Glial Cells: Recent Advances. Springer Singapore: Singapore; 2022. p. 401–37.

    Chapter  Google Scholar 

  92. Kang HJ, Kawasawa YI, Cheng F, Zhu Y, Xu X, Li M, et al. Spatio-temporal transcriptome of the human brain. Nature. 2011;478(7370):483–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Liu L, Zhang Y, Niu G, Li Q, Li Z, Zhu T, et al. BrainBase: a curated knowledgebase for brain diseases. Nucleic Acids Res. 2021;50(D1):D1131–8.

    Article  PubMed Central  Google Scholar 

  94. Hochgerner H, Singh S, Tibi M, Lin Z, Skarbianskis N, Admati I, et al. Neuronal types in the mouse amygdala and their transcriptional response to fear conditioning. Nat Neurosci. 2023;26(12):2237–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Kunimatsu J, Maeda K, Hikosaka O. The Caudal Part of Putamen Represents the Historical Object Value Information. J Neurosci. 2019;39(9):1709–19.

    CAS  PubMed  PubMed Central  Google Scholar 

  96. Gordji-Nejad A, Matusch A, Li S, Kroll T, Beer S, Elmenhorst D, et al. Phosphocreatine levels in the left thalamus decline during wakefulness and increase after a nap. J Neurosci. 2018;38(49):10552.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  97. Hayasaka N, Aoki K, Kinoshita S, Yamaguchi S, Wakefield JK, Tsuji-Kawahara S, et al. Attenuated food anticipatory activity and abnormal circadian locomotor rhythms in Rgs16 knockdown mice. PLoS ONE. 2011;6(3): e17655.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Lipiec MA, Bem J, Koziński K, Chakraborty C, Urban-Ciećko J, Zajkowski T, et al. TCF7L2 regulates postmitotic differentiation programmes and excitability patterns in the thalamus. Development. 2020;147(16):dev190181.

  99. Chen D-H, Méneret A, Friedman JR, Korvatska O, Gad A, Bonkowski ES, et al. ADCY5-related dyskinesia. Neurology. 2015;85(23):2026–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Waheed A, Sly WS. Carbonic anhydrase XII functions in health and disease. Gene. 2017;623:33–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Hua K, Wang P, Lan Z, Li M, Zhao W, Wang T, et al. Increased left putamen volume correlates with pain in ankylosing spondylitis patients. Front Neurol. 2020;11: 607646.

    Article  PubMed  PubMed Central  Google Scholar 

  102. Blank MC, Grinberg I, Aryee E, Laliberte C, Chizhikov VV, Henkelman RM, et al. Multiple developmental programs are altered by loss of Zic1 and Zic4 to cause Dandy-Walker malformation cerebellar pathogenesis. Development. 2011;138(6):1207–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Juan SMA, Adlard PA. Ageing and cognition. In: Harris JR, Korolchuk VI, editors. Biochemistry and Cell Biology of Ageing: Part II Clinical Science. Springer Singapore: Singapore; 2019. p. 107–22.

    Chapter  Google Scholar 

  104. Liu J, Zhang B, Lei H, Feng Z, Liu J, Hsu A-L, et al. Functional aging in the nervous system contributes to age-dependent motor activity decline in C. elegans. Cell Metab. 2013;18(3):392–402.

  105. Ximerakis M, Lipnick SL, Innes BT, Simmons SK, Adiconis X, Dionne D, et al. Single-cell transcriptomic profiling of the aging mouse brain. Nat Neurosci. 2019;22(10):1696–708.

    Article  CAS  PubMed  Google Scholar 

  106. Hajdarovic KH, Yu D, Hassell L-A, Evans SA, Packer S, Neretti N, et al. Single-cell analysis of the aging female mouse hypothalamus. Nat Aging. 2022;2(7):662–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Zhang W, Zhang S, Yan P, Ren J, Song M, Li J, et al. A single-cell transcriptomic landscape of primate arterial aging. Nat Commun. 2020;11(1):2202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  108. Iskusnykh IY, Zakharova AA, Kryl'skii ED, Popova TN. Aging, Neurodegenerative Disorders, and Cerebellum. Int J Mol Sci. 2024;25(2).

  109. Campeau A, Mills RH, Stevens T, Rossitto L-A, Meehan M, Dorrestein P, et al. Multi-omics of human plasma reveals molecular features of dysregulated inflammation and accelerated aging in schizophrenia. Mol Psychiatry. 2022;27(2):1217–25.

    Article  CAS  PubMed  Google Scholar 

  110. Eavri R, Shepherd J, Welsh CA, Flanders GH, Bear MF, Nedivi E. Interneuron Simplification and Loss of Structural Plasticity As Markers of Aging-Related Functional Decline. J Neurosci. 2018;38(39):8421–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Hanslik KL, Marino KM, Ulland TK. Modulation of Glial Function in Health, Aging, and Neurodegenerative Disease. Front Cell Neurosci. 2021;15.

  112. Atlas A. a multi-omics database for aging biology. Nucleic Acids Res. 2021;49(D1):D825-d830.

    Article  Google Scholar 

  113. Jia K, Cui C, Gao Y, Zhou Y, Cui Q. An analysis of aging-related genes derived from the Genotype-Tissue Expression project (GTEx). Cell Death Discov. 2018;4(1):91.

    Article  PubMed  PubMed Central  Google Scholar 

  114. Jin LE, Wang M, Yang ST, Yang Y, Galvin VC, Lightbourne TC, et al. mGluR2/3 mechanisms in primate dorsolateral prefrontal cortex: evidence for both presynaptic and postsynaptic actions. Mol Psychiatry. 2017;22(11):1615–25.

    Article  CAS  PubMed  Google Scholar 

  115. Guo Z, Zhang J, Liu X, Hou H, Cao Y, Wei F, et al. Neurometabolic characteristics in the anterior cingulate gyrus of Alzheimer’s disease patients with depression: a 1H magnetic resonance spectroscopy study. BMC Psychiatry. 2015;15(1):306.

    Article  PubMed  PubMed Central  Google Scholar 

  116. Kant S, Stopa EG, Johanson CE, Baird A, Silverberg GD. Choroid plexus genes for CSF production and brain homeostasis are altered in Alzheimer’s disease. Fluids Barriers CNS. 2018;15(1):34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  117. Fujikawa M, Ohsakaya S, Sugawara K, Yoshida M. Population of ATP synthase molecules in mitochondria is limited by available 6.8-kDa proteolipid protein (MLQ). Genes Cells. 2014;19(2):153–160.

  118. Zuehlke AD, Beebe K, Neckers L, Prince T. Regulation and function of the human HSP90AA1 gene. Gene. 2015;570(1):8–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  119. Hu B, Zhang Y, Jia L, Wu H, Fan C, Sun Y, et al. Binding of the pathogen receptor HSP90AA1 to avibirnavirus VP2 induces autophagy by inactivating the AKT-MTOR pathway. Autophagy. 2015;11(3):503–15.

    Article  PubMed  PubMed Central  Google Scholar 

  120. Braunewell KH, Klein-Szanto AJ. Visinin-like proteins (VSNLs): interaction partners and emerging functions in signal transduction of a subfamily of neuronal Ca2+ -sensor proteins. Cell Tissue Res. 2009;335(2):301–16.

    Article  CAS  PubMed  Google Scholar 

  121. Lin CW, Chang LC, Tseng GC, Kirkwood CM, Sibille EL, Sweet RA. VSNL1 Co-Expression Networks in Aging Include Calcium Signaling, Synaptic Plasticity, and Alzheimer’s Disease Pathways. Front Psychiatry. 2015;6:30.

    Article  PubMed  PubMed Central  Google Scholar 

  122. Braunewell KH, Dwary AD, Richter F, Trappe K, Zhao C, Giegling I, et al. Association of VSNL1 with schizophrenia, frontal cortical function, and biological significance for its gene product as a modulator of cAMP levels and neuronal morphology. Transl Psychiatry. 2011;1(7):e22–e22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  123. Li H, Dai Y, Luo Z, Nie D. Cloning of a new testis-enriched gene C4orf22 and its role in cell cycle and apoptosis in mouse spermatogenic cells. Mol Biol Rep. 2019;46(2):2029–38.

    Article  CAS  PubMed  Google Scholar 

  124. Damaj L, Lupien-Meilleur A, Lortie A, Riou É, Ospina LH, Gagnon L, et al. CACNA1A haploinsufficiency causes cognitive impairment, autism and epileptic encephalopathy with mild cerebellar symptoms. Eur J Hum Genet. 2015;23(11):1505–12.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  125. de Oliveira Figueiredo EC, Calì C, Petrelli F, Bezzi P. Emerging evidence for astrocyte dysfunction in schizophrenia. Glia. 2022;70(9):1585–604.

    Article  PubMed  PubMed Central  Google Scholar 

  126. Ahrendts J, Rüsch N, Wilke M, Philipsen A, Eickhoff SB, Glauche V, et al. Visual cortex abnormalities in adults with ADHD: a structural MRI study. World J Biol Psychiatry. 2011;12(4):260–70.

    Article  PubMed  Google Scholar 

  127. Yokokura M, Takebasashi K, Takao A, Nakaizumi K, Yoshikawa E, Futatsubashi M, et al. In vivo imaging of dopamine D1 receptor and activated microglia in attention-deficit/hyperactivity disorder: a positron emission tomography study. Mol Psychiatry. 2021;26(9):4958–67.

    Article  CAS  PubMed  Google Scholar 

  128. Nagayasu Y, Ito J-i, Nishida T, Yokoyama S. Reactivity of astrocytes to fibroblast growth factor-1 for biogenesis of apolipoprotein E-high density lipoprotein is down-regulated by long-time secondary culture. J Biochem. 2008;143(5):611–616.

  129. Park H-J, Friston K. Structural and functional brain networks: From connections to cognition. Science. 2013;342(6158):1238411.

    Article  PubMed  Google Scholar 

  130. Wang Y-M, Wang W-C, Pan Y, Zeng L, Wu J, Wang Z-B, et al. Nucleus versus gene UMI count matrix across multiple distinct non-human primate brain regions from single-nucleus RNA sequencing. Science Data Bank. https://doiorg.publicaciones.saludcastillayleon.es/10.57760/sciencedb.06819. 2022.

  131. Chen T, Chen X, Zhang S, Zhu J, Tang B, Wang A, et al. The Genome Sequence Archive family: Toward explosive data growth and diverse data types. Genom Proteom Bioinf. 2021;19(4):578–83.

    Article  Google Scholar 

  132. Wang Y-M, Wang W-C, Pan Y, Zeng L, Wu J, Wang Z-B, et al. Cellular architecture of non-human primate brains. Dataset CRA009154. Genome Sequence Archive. https://bigd.big.ac.cn/gsa/browse/CRA009154. 2025.

  133. Wang Y-M, Wang W-C, Pan Y, Zeng L, Wu J, Wang Z-B, et al. Regional and aging-specific cellular architecture of non-human primate brains. GitHub. https://github.com/WanYMEN/Regional-and-aging-specific-cellular-architecture-of-non-human-primate-brains. 2025.

Download references

Acknowledgements

We would like to thank the Institutional Center for Shared Technologies and Facilities of Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences (CAS) for providing us with the confocal microscope. We are also grateful to Guolan Ma for her technical support and Zhengfei Hu of KIZ, CAS, for his assistance with animal care.

Funding

This study was funded by the National Key R&D Program of China (2023YFA1800500, 2021YFF0702703), the STI2030-Major Projects (2021ZD0200900), the Key-Area Research and Development Program of Guangdong Province (2019B030335001), Yunnan Fundamental Research Projects (202305AH340006), the National Key R&D Program of China (2018YFA0801403, 2021YFF0702700), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB32060200), the National Natural Science Foundation of China (81941014, 81960422), the STI2030-Major Projects (2022ZD0212700), the Yunnan Fundamental Research Projects (202201AT070139, 202201AT070153), the Applied Basic Research Programs of Science Technology Commission Foundation of Yunnan Province (202001AT070130), and the Science and Technology Project of Yunnan Province (202101AY070001- 001).

Author information

Authors and Affiliations

Authors

Contributions

D.-D.W., X.-T.H., and Y.H. conceived and supervised the study. Y.-M.W. performed the data analysis. W.-C.W., J.W., L.-M.W., and Y.-Y.F. performed the immunofluorescence and immunohistochemistry staining experiments. Y.P. conducted all experiments associated with knockdown of VSNL1 and HPCAL4 in U251 cells. Y.-M.W., L.Z., and Z.-B.W. collected the samples. X.-L.Z. prepared the samples for snRNA-seq. Z.-B.W. measured the spatial distance between regions. M.-L.L., S.W., and Y.S. provided useful assistance with data analysis. D.-D.W., Y.-M.W., W.-C.W., and Y.P. wrote the manuscript. X.-T.H., Y.H., and D.N.C. revised the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yonghan He, Xin-Tian Hu or Dong-Dong Wu.

Ethics declarations

Ethics approval and consent to participate

All experimental protocols followed the animal experimentation guidelines and regulations of the Kunming Institute of Zoology, Chinese Academy of Sciences (Yunnan, China). This research was reviewed and approved by the Institutional Animal Care and Use Committee of the Kunming Institute of Zoology (permit nos. IACUC-PE-2022-01- 003).

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

13073_2025_1469_MOESM1_ESM.pdf

Additional file 1: Fig. S1. Anatomy of the ten dissected brain regions used for snRNA-seq in rhesus macaque. Fig. S2. Statistical information summarizing the data across the 39 brain samples. Fig. S3. Characterization of consistent expression profiles. Fig. S4. Comparisons of our cell annotation with published NHP brain snRNA-seq datasets. Fig. S5. Subtype features of the putamen-specific group. Fig. S6. Subtype features of the cerebellum-specific group. Fig. S7. Subtype features of the thalamus-specific group. Fig. S8. Subtype features of the excitatory group. Fig. S9. Subtype features of the inhibitory group. Fig. S10. Subtype features of the astrocytic group. Fig. S11. Subtype features of the microglial group. Fig. S12. Subtype features of the oligodendrocyte-lineage group. Fig. S13. Subtype features of the membrane-related non-neuronal cell type group. Fig. S14. Validation of protein expression across the ten brain regions of interest for the selected region-specific genes. Fig. S15. Examination of data distributions of the selected RSGs from immunofluorescence staining. Fig. S16. Transcriptomic features across brain regions. Fig. S17. Correlations between intra- or inter-region cellular communication and spatial distance. Fig. S18. Boxplot indicates comparison of transcriptional noise across different cell types and brain regions between young and old individuals. Fig. S19. Alterations in cellular communication with aging. Fig. S20. Comparison of gene expression profiles between young and old macaque brains. Fig. S21. Examination of data distributions of the selected aging-related genes from immunohistochemical staining. Fig. S22. Expression of HPCAL4 and VSNL1 across cell types. Fig. S23. Knockdown of VSNL1 and HPCAL4 in U251 glioma cells. Fig. S24. Examination of data distributions after knockdown of VSNL1 or HPCAL4 in U251 glioma cells. Fig. S25. Inferred progressive aging trajectories across cell types or subtypes using monocle3. Fig. S26. Heritability enrichment of human traits/diseases across signatures related to regions and aging.

Additional file 2: Table S1. Information of samples used for snRNA-seq and protein validation.

Additional file 3: Table S2. List of cell type and subtype marker genes.

Additional file 4: Table S3. List of region-specific genes.

Additional file 5: Table S4. List of DEGs between young and old macaques across cell types and regions.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, YM., Wang, WC., Pan, Y. et al. Regional and aging-specific cellular architecture of non-human primate brains. Genome Med 17, 41 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13073-025-01469-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13073-025-01469-x

Keywords