Lifespan single-cell transcriptomic atlas of the human prefrontal cortex

DLPFC lifespan study design

Brain tissue specimens were obtained from NIMH-IRP Human Brain Collection Core (HBCC) (172 samples, from individuals aged 0.2–85 years) and The Mount Sinai NIH Neurobiobank (MSSM) (112 samples, from individuals aged 20–97 years). In total, 284 neurotypical controls aged 0–97 years from the PsychAD dataset were included in this study (Supplementary Fig. 1a). The majority of the samples are from individuals of European (n = 158) followed by African (n = 95), American (n = 26), East Asian and South Asian (n = 5) descent. Supplementary Fig. 1a shows the demographic information at the donor level, including sex, age, time of death and ancestry, stratified by corresponding brain banks. The data were categorized into four groups: (1) developmental, which constitutes infancy (0–1 year), childhood (2–11 years) and adolescence (12–19 years); and (2) young (20–39 years), (3) middle (40–59 years) and (4) late adulthood (≥60 years).

We used the available neuropathology details on MSSM samples for control sample selection. The selection criterion for MSSM included the following:

1. (1)

   CERAD scores: for neuritic plaque density of MSSM samples, only those with a score of 1 (no AD) were included.

2. (2)

   Braak stage: samples with Braak stages 0, 1 or 2 were retained.

3. (3)

   Secondary diagnoses: donors with any additional brain-related diagnosis, including neurodegenerative (for example, AD and Parkinson’s disease) and neuropsychiatric diseases (for example, schizophrenia and bipolar disorder) as well as the presence of mild cognitive impairment, were not retained.

In principle, we applied equivalent selection criteria to the HBCC samples. Although the HBCC cohort did not provide specific Braak and CERAD values, we confirmed through review of neuropathological reports that the selected donors did not exhibit substantial plaque or tangle pathology. Thus, all selected donors who lacked brain-related diagnoses were deemed reliable neurotypical controls, despite the absence of detailed neuropathological data.

The link to the complete demographic and clinical information of the present study population is provided in Supplementary Table 1 and Supplementary Data 1.

FANS protocol and snRNA-seq hashing from frozen brain tissue

All libraries from the PsychAD dataset were prepared using a standardized protocol for nuclei isolation and hashing¹⁸. The dataset was curated after completing snRNA-seq preprocessing and the taxonomy step, as described in this section and the following sections. Library generation involved isolating and sorting nuclei from frozen brain specimens using fluorescence-activated nuclear sorting (FANS). Then, 25 mg of frozen post-mortem human brain tissue was homogenized in a cold lysis buffer with RNase inhibitors. The homogenate was filtered through a 40 µm cell strainer, and the flow-through was underlaid with sucrose solution before centrifugation at 107,000g for 1 h at 4 °C. The resulting pellets were resuspended in PBS with 0.5% BSA. Six samples were processed simultaneously, with up to 2 million nuclei per sample pelleted at 500g for 5 min at 4˚ C. The nuclei were then resuspended in 100 µl staining buffer and incubated with 1 µg of a unique TotalSeq-A nuclear hashing antibody (BioLegend) for 30 min at 4 °C. Before FANS, the volumes were adjusted to 250 µl with PBS, and 7-aminoactinomycin D (7-AAD) was added according to the manufacturer’s instructions. The 7-AAD-positive nuclei were sorted into tubes precoated with 5% BSA using the FACSAria flow cytometer (BD Biosciences). FACSDiva software (BD Biosciences, v.8.0.2) was used for data collection.

After FANS, the nuclei were washed twice with 200 µl of staining buffer, resuspended in PBS and quantified using the Countess II (Life Technologies). The concentrations were adjusted, and equal volumes of differentially hash-tagged nuclei were combined. Using 10x Genomics single cell 3′ v3.1 reagents, 60,000 nuclei (10,000 per donor) were processed in each of two lanes to create technical replicates. During cDNA amplification for library preparation, 1 µl of a 2 µm HTO cDNA PCR additive primer⁵⁰ was included. The supernatant from a 0.6× SPRI selection was reserved for HTO library generation. Both cDNA and HTO libraries were prepared according to the manufacturer’s instructions (10x Genomics and BioLegend, respectively). Sequencing was performed at the NYGC using the NovaSeq platform (Illumina).

Preprocessing of the snRNA-seq dataset

The DLPFC lifespan cohort analysed here is a subset of the broader PsychAD dataset^8,9. Below, we summarize the computational processing and quality control applied to the full PsychAD dataset, followed by the dimensionality reduction steps that were recomputed after restricting analyses to the DLPFC lifespan cohort.

Computational processing and demultiplexing

Sequencing reads from all pools of multiplexed samples were aligned to the hg38 reference genome using STARsolo (v.2.7.9a)⁵¹. To assign cells from each pool to their respective donors, a genotype-based demultiplexing pipeline was employed, followed by genotype concordance checks. Specifically:

1. (1)

   Pile-up of alleles: using cellSNP (v.1.2.0)⁵², alleles were aggregated from polymorphic sites overlapping snRNA-seq reads within expressed genes (genes expressed by at least 10 cells were included). Polymorphic sites required a minimum minor allele frequency of 0.1 and a minimum aggregated unique molecular identifier (UMI) count of 20.

2. (2)

   Donor assignment: Vireo (v.0.5.8)⁵³ was used to cluster cells into six distinct donor groups per pool based on allele pile-ups. Each cluster’s identity was assigned through genotype concordance analysis, comparing cell clusters to reference genotyping data using QTLtools-mbv (v.1.3)⁵⁴.

3. (3)

   Baseline quality control for genotyping: only cells meeting baseline quality-control thresholds (expressing at least 1,000 genes and with a mitochondrial read fraction below 5%) were included in this analysis.

This pipeline detected and corrected occasional sample swaps and mislabelling, ensuring accurate donor assignment for the majority of pools. These steps were performed once on the full PsychAD dataset^8,9 before any downstream subsetting.

Rigorous quality-control pipeline

After genome alignment and demultiplexing, downstream quality control and processing were conducted using Pegasus (v.1.7.0)⁵⁵ and Scanpy (v.1.9.1)⁵⁶, with additional steps to ensure robust data integrity:

1. (1)

   Individual cell quality control: cells were excluded if their UMI counts were outside the range of 1,500 ≤ UMI ≤ 110,000, gene counts were outside 1,100 ≤ genes ≤ 12,500, or if their mitochondrial content exceeded 5% (mitochondrial fraction > 0.05). Ambient RNA contamination was mitigated using CellBender (v.0.2.1)⁵⁷, as well as the proportion of reads mapped to non-mRNA categories like rRNA, sRNA and pseudogenes, in addition to examining confounding factors such as the lncRNA MALAT1. Doublets were identified and removed using Scrublet (v.0.2.3)⁵⁸.

2. (2)

   Feature-level quality control: genes not expressed in at least 0.05% of nuclei were removed.

3. (3)

   Donor-level quality control: donors with fewer than 50 nuclei were excluded to avoid introducing noise into downstream analyses.

These quality-control filters were applied uniformly across the full PsychAD dataset^8,9 and were not reapplied after subsetting to the DLPFC lifespan cohort.

Data integration

To correct for (non-biological) variance in the full PsychAD dataset^8,9, including tissue dissection biases or brain bank resources, we performed canonical correlation analysis using Harmony (v.0.1)⁵⁹. After identifying highly variable features, a k-nearest-neighbour graph was constructed using Harmony-corrected principal component analysis (PCA) embeddings, and the Leiden algorithm was applied to cluster cells by type. Clustering annotations were retained for downstream analyses unless otherwise stated.

Basic quality-control metrics of the DLPFC lifespan dataset are presented in Supplementary Fig. 2a. These include the distribution of nuclei across donors and a histogram of nuclei per sample, as well as distribution plots for mitotic ratio, ribosomal gene content, total counts (n_counts) and the number of genes (n_genes), all split by brain bank.

Joint cellular taxonomy

Cellular taxonomy for the DLPFC lifespan cohort was retained from the PsychAD dataset^8,9, which includes comprehensive classification at the class, subclass and subtype levels. Cell-type labels for the DLPFC lifespan cohort were inherited directly from the PsychAD reference atlas and were not recomputed. The PsychAD taxonomy was established using a modular and iterative clustering approach^8,9 with marker validation and label transfer using scANVI (v.0.20.3)⁶⁰. The validation of taxonomy was done using the reference dataset from ref. ¹⁰ and through Xenium spatial experiments. The hierarchy comprises class, subclass and subtype derived from a previous study¹⁰. Supplementary Fig. 3 presents the subclass-level taxonomy pseudobulk correlation between DLPFC lifespan nuclei and the reference dataset¹⁰, as well as the expression of representative marker genes.

UMAP

To ensure an appropriate low-dimensional representation after restricting analyses to the DLPFC lifespan cohort, dimensionality reduction was recomputed. To reduce memory usage and standardize the feature space before recomputation, an on-disk subset of the .h5ad file was generated using custom ondisk_subset function (chunk_size = 500,000), retaining all nuclei and restricting genes to protein-coding autosomal features (gene_type == “protein_coding” and excluding chromosomes MT, X and Y); the .raw slot was preserved (raw = True).

Highly variable genes (HVGs) were then identified from the subsetted .h5ad using a Scanpy-based HVG procedure implemented in custom function scanpy_hvf_h5ad with flavor = “cell_ranger” and batch-aware selection (batch_key = “Source”). Mean–dispersion thresholds were applied (min_mean = 0.0125, max_mean = 3, min_disp = 0.5), and the top 6,000 HVGs were retained (n_top_genes = 6000). This HVG list was imported into Pegasus⁵⁵ (v.1.10.1) using the setting data.var[“highly_variable_features”].

Dimensionality reduction was performed in Pegasus using PCA (pg.pca) with up to 30 components. Before batch correction, total UMI counts (n_counts), mitochondrial read fraction (percent_mito) and cell cycle effects (cycle_diff) were regressed using pg.regress_out. Batch correction was then applied using Harmony via pg.run_harmony (batch = “Source”, rep = “pca_regressed”, max_iter_harmony = 20), retaining 30 components. A k-nearest-neighbour graph was constructed using pg.neighbours on the Harmony-corrected PCA embedding (rep = “pca_regressed_harmony”, K = 100, dist = “l2”, n_comps = 30). UMAP embedding was subsequently computed using pg.umap on the Harmony-corrected embedding (rep = “pca_regressed_harmony”, n_neighbors = 100, rep_ncomps = 30) (Fig. 1a). UMAT calculation (Supplementary Fig. 1b) was run based on the approach introduced previously⁶.

Single-cell polygenic disease risk score

We used the scDRS package (v.1.0.1)⁵⁶ to evaluate the combined expression of potential disease-associated genes obtained from GWAS summary statistics through MAGMA analysis⁶¹ (Supplementary Fig. 4). Each gene’s contribution was weighted by its MAGMA z score from GWAS and inversely weighted by its gene-specific technical noise level in the single-cell data. This process was performed across each cell of the DLPFC lifespan dataset to generate raw disease scores specific to each nucleus (scdrs.preprocess; n_mean_bin = 20, n_var_bin = 20). Moreover, we generated 200 sets of raw control scores, matched in gene set size, mean expression and expression variance to the disease-associated genes. Subsequently, we normalized both the raw disease scores and raw control scores for each cell, resulting in normalized scores. These calculations were executed using the scdrs.score_cell function, with the following parameters: scdrs.score_cell(–ctrl_match_key = “mean_var”,–n_ctrl=200,–weight_opt = “vs”,–return_ctrl_raw_score=False,–return_ctrl_norm_score=True,–verbose=False).

For further analysis, subclass-level examinations were conducted to link predefined subclasses to disease and evaluate the heterogeneity in disease association across cells within each predefined subclass level of taxonomy. This was achieved using the scdrs.method.downstream_group_analysis function with the default settings. The output from this step is provided in Supplementary Data 2. A list of all 24 GWAS traits used in the study is provided in Supplementary Table 2.

Pseudobulk data aggregation and covariate selection

To quantify the variance explained by age, age-associated changes in the transcriptome and lifespan trends, we pseudobulked gene expression data by aggregating the counts for each donor from 1,307,674 nuclei, stratified by subclass and four age groups (developmental, young, middle and late adulthood) using the aggregateToPseudoBulk function from the dreamlet (v.1.1.17) R package. Subsequently, we applied voom normalization to the age-group-subclass-specific gene-by-donor matrix using the processAssays function, which filters for samples with at least 5 nuclei per subclass and genes with a minimum of 5 reads per donor. This process generated four lists containing gene × donor matrices for each subclass corresponding to each age group.

With source, sex and PMI (post-mortem interval in hours) as the base model, technical and biological covariates were identified on the entire PsychAD dataset. As the DLPFC lifespan nuclei are a subset of the PsychAD dataset, we implemented the same model for i subclass and the j gene in our analysis as described in equation (1). Hereafter, we use these variables as covariates.

The details of these covariates are as follows: PMI is coded as a scaled numerical variable, while source and sex are categorical variables. Source is excluded from the equation for any downstream analyses for the developmental group as all samples are from a single brain bank (HBCC) (Supplementary Data 1). The number of genes from each nucleus, n genes, and the proportion of mitochondrial genes per donor, percentage mito, are summarized for each donor and were obtained using the qc_metrics function from the pegasus (v.1.8.1)⁵⁵ package. Scores for reference mitochondrial genes from the mitochondrial chromosome, mito genes, and mitochondrial ribosomal genes not from the mitochondrial chromosome, mito ribo, as well as ribosomal genes from both large and small units, ribo genes, are estimated for each nucleus and summarized for each donor using the calc_signature_score function from pegasus (v.1.8.1)⁵⁵. Supplementary Fig. 2c shows the pairwise correlation of these covariates with age.

Quantification of variance in age

We quantified the variance explained by age for age-group-subclass-specific pseudobulk expression data for i subclass and gene j as expression_ij = age + covariates in the fitVarPart function from the dreamlet (v.1.1.17) R package. VariancePartition¹¹ is a statistical framework designed to prioritize drivers of variation in gene expression, using a linear mixed model to quantify contributions from factors such as disease status, sex, cell type and technical variables. The covariates were obtained from equation (1). Fig. 1c highlights the mean contribution of numerical variable age in pseudobulk expression for each subclass across four age groups. Supplementary Fig. 2d shows the combined contribution of covariates and age from four groups. The output from this section is provided in Supplementary Data 3.

Lifespan dynamics of nucleus counts

To quantify changes in nucleus counts as a function of age, we used the crumblr function from the crumblr¹⁵ (count ratio uncertainty modelling based linear regression) (v.0.99.6) R package. Traditional methods (such as per-cell-type proportion tests) may miss subtle, coordinated shifts across related subtypes and often overlook covariate effects. crumblr addresses this by applying precision-weighted linear mixed models within a multivariate framework, improving sensitivity to ageing-related cell type composition changes across the cell lineage hierarchy. The package uses a three-step process to quantify the association of changes in nucleus counts with user-defined independent variables. (1) Normalization: nucleus counts are normalized using a Dirichlet multinomial distribution. (2) Modelling: a standard dream precision-weighted linear mixed model with empirical Bayes estimation is implemented to obtain association statistics for each measurement. (3) Multivariate hypothesis testing: this step enables the joint analysis of internal nodes in a hierarchical clustering of subclasses, which is crucial for accounting for correlations among related subclasses, such as those among EN subclasses. Sample code for these steps is provided in the repository (Code availability).

In this section, we performed two analyses to quantify changes in nucleus counts: (1) as a function of age across the entire lifespan; and (2) within specific age groups. These analyses were limited to subclasses with at least 500 nuclei counts resulting in 26 out of 27 subclasses. For both analyses, we explored the optimal model to represent nucleus counts as a function of age, aiming to determine whether a linear or logarithmic relationship would be optimal. To test these models, we first regressed out the subset of covariates: sex, PMI and source from nucleus counts for i subclass across 284 donors using the equations (2) and (3). The other covariates are more relevant to gene expression analysis; therefore, we excluded them from the nucleus composition analysis.

After this, we tested two models:

Using difference of BIC (ΔBIC_i = BIC₂ − BIC₁), for the i subclass, we found that log₂[age + 1] showed an optimal relationship between age and changes in nucleus counts for most subclasses as shown in Supplementary Fig. 6a. Applying k-means clustering to subclass-specific coefficients from the log₂[age + 1] model (Supplementary Fig. 6b), we identified two distinct clusters: one with positive coefficients indicating a log-increasing trend, and another with negative coefficients indicating a log-decreasing trend, as illustrated in Extended Data Fig. 1a. To further quantify these lifespan associations, we performed crumblr-enabled modelling and multivariate test analysis using equation (4). The coefficient of log₂[age + 1] is shown in Extended Data Fig. 1c and the output is provided in Supplementary Table 3.

Next, we performed age-group-specific analysis by first quantifying the variance explained by age in normalized nucleus counts for each age group using fitExtractVarPartModel from the variancePartition¹¹ (v.1.33.11) R package including equation (4) as a model. Figure 1e shows the distribution of variance explained by log₂[age + 1] in 26 subclasses from four age groups. Subsequently, we quantified the coefficients of log₂[age + 1] and performed crumblr-enabled modelling and multivariate hypothesis test using equation (4) for each age group (the source covariate was removed from equation (4) for the developmental group). Figure 1f shows the heat map of the coefficient of log₂[age + 1] for each subclass and age group. The results of this analysis are provided in Supplementary Table 3.

To test whether early developmental changes (0–20 years) may dominate the modelling, we performed a breakpoint analysis (Supplementary Fig. 5c,d). We iteratively split the dataset at age thresholds ranging from 20 to 30 years and modelled changes in nucleus counts separately for the

Finally, to investigate potential biological mechanisms underlying the most vulnerable age group to compositional changes in nuclei, we performed a pathway activity analysis stratified by subclasses showing log-increasing or log-decreasing age trends. We first curated a set of 374 brain-relevant GO biological process pathways from MSigDB using keyword-based filtering with terms such as neuron, glia, synapse, axon, oligodendrocyte, astrocyte, microglia and cortex. For each nucleus, pathway activity scores were calculated using the over-representation analysis method implemented in the decoupler (v.1.8.0)³¹ Python package. This method performs over-representation analysis using a Fisher’s exact test applied to the top 5% most highly expressed genes in each nucleus, relative to the genes annotated to each pathway. The resulting nucleus-by-pathway activity score matrix was aggregated at the subclass level by averaging across all nuclei within each subclass. Subclasses were then annotated with their lifespan trajectory group (log-increasing or log-decreasing), as defined by the compositional modelling described above. To identify biological pathways associated with these ageing trends, we performed Welch’s t-tests comparing subclass-level activity scores between the two trajectory groups. Pathways with FDR-adjusted P < 0.05 were retained and ranked by the absolute difference in mean activity scores. We selected the top 15 pathways with higher activity in log-increasing subclasses and the top 15 with higher activity in log-decreasing subclasses for visualization. Pathway scores were min–max normalized across subclasses, and subclass annotations were used to group columns by major cell class as shown in the heat map (Extended Data Fig. 1f)

Validation of nucleus composition changes with age using RNAscope

To validate age-associated changes in nucleus composition, we collected formalin-fixed paraffin-embedded (FFPE) tissue blocks from the DLPFC of eight donors across the lifespan. Samples were obtained from the Neuropathology Brain Bank and Research CoRE at the Icahn School of Medicine at Mount Sinai. Donors were stratified into four life stages: developmental (n = 3): 0.33, 0.75 and 9 years old; young adulthood (n = 1): 38 years old; middle adulthood (n = 2): 51 and 53 years old; and late adulthood (n = 2): 65 and 86 years old. To assess subclass-specific changes in nucleus composition across the lifespan, we selected marker genes for the IN_SST and oligo subclasses, which exhibited opposing age-related trends: log-linear decreases for IN_SST and increases for oligo (Extended Data Fig. 1e). IN_SST cells were labelled using the SST gene, detected using the RNAscope probe Hs-SST-C2 in channel C2 (310591-C2), while oligodendrocytes were labelled using the GPR37 gene, detected using the RNAscope probe Hs-GPR37-C3 in channel C3 (513631-C3).

RNAscope assay and TrueBlack treatment pipeline

RNAscope Multiplex Fluorescent v2 assay (ACD, UM 323100/Rev B) was performed on FFPE DLPFC sections (5 µm). Slides were deparaffinized, rehydrated and underwent target retrieval using preheated 1× target retrieval reagent (≥99 °C), followed by ethanol dehydration and drying. A hydrophobic barrier was drawn around tissue sections, which were then treated with hydrogen peroxide (10 min) and protease IV (30 min) at room temperature. Gene-specific probes were hybridized at 40 °C for 2 h, followed by AMP1–3 incubations (30 min each) with washes in between. Fluorescence signal amplification was performed using sequential HRP-C1, C2, and C3 detection, each with distinct fluorophores. Nuclei were counterstained with DAPI. To reduce lipofuscin autofluorescence, slides were treated with 1× TrueBlack (23007, Biotium) in 70% ethanol for 30 s, then washed with PBS and mounted using ProLong Gold Antifade Mountant (P36930, Thermo Fisher Scientific).

Imaging

RNAscope images were acquired using the EVOS M7000 microscope (Thermo Fisher Scientific) using a ×40 objective (OLY XAPO ×40/0.95 NA). For each tissue section, 3–4 high-quality regions of interest (ROIs) were selected to capture representative areas spanning the full grey matter and half of the white matter. Imaging was performed using the DAPI, RFP and Cy5 filter cubes to detect nuclei and RNAscope signals (C2 and C3 channels). ROIs were scanned using the serpentine horizontal pattern in quick scan mode with autofocus on every DAPI field. Each ROI consisted of 100–200 fields of view (FOVs), saved as 16-bit TIFF files.

Image analysis

FOVs were stitched into full ROI images using a custom Fiji (ImageJ) macro based on the Grid/Collection Stitching Plugin. DAPI channels were stitched first (snake-by-rows layout with subpixel accuracy), generating a reference TileConfiguration.registered.txt used to stitch the corresponding C2 and C3 channels for alignment. Stitched channel images were combined into hyperstacks using Fiji’s Images to Stack and Stack to Hyperstack functions. For quantification, stitched ROIs were analysed in QuPath. Cells were segmented using the Cell Detection tool, and RNAscope puncta were detected with the Subcellular Detection function. Single-cell expression and composition data were exported using saveDetectionMeasurements.

Statistical analysis

To quantify subclass-specific expression, we analysed QuPath output tables from stitched RNAscope images across eight donors. For each ROI, we calculated the fraction of positive cells by dividing the number of cells with ≥2 RNA molecules for the target marker (for example, SST) by the total number of detected nuclei in that ROI. This metric quantifies the relative abundance of marker-positive cells and serves as a proxy for cell composition across age groups. Similar to the lifespan model as described in Extended Data Fig. 1a, the age of eight donors was modelled as log₂[age + 1]. Statistical analysis and visualizations were performed using the dplyr (v.1.1.4), ggplot2 (v.4.0.0) R packages.

Age-associated transcriptomic changes

To assess the age-related transcriptomic changes for each subclass and age group, we performed dreamlet enabled limma modelling on age group-subclass-specific pseudobulk expression data as described in the ‘Pseudobulk data aggregation and covariate selection’ section. The model implemented for i subclass and the j gene is expression_ij = age + covariates. We conducted final multiple testing correction at the study-wide level (26 subclasses × number of genes per subclass) for each group. The total number of study-wide genes from all subclasses was as follows: development (n = 328,991), young adulthood (n = 333,471), middle adulthood (n = 340,084) and late adulthood (n = 292,748). After applying a threshold of FDR < 0.05, we obtained aDEGs per age group, as shown in Fig. 1h and Supplementary Fig. 7b,c. We provide the results of this analysis in Supplementary Data 4. For gene set pathway analysis to identify biological processes, we used the enrichr function from the enrichR (v.3.2) R library. The database used was GO Biological Process 2021, and the results are shown in Fig. 1i and Supplementary Fig. 7f. The pathway enrichment analysis for all aDEGs per subclass is provided in Supplementary Data 5. dreamlet enabled us to model donor-level variation and batch effects using precision-weighted mixed models, outperforming standard methods. Its efficiency and design-aware approach ensured statistically rigorous identification of DEGs.

Lifespan trends of gene expression

To obtain the lifespan trend of each gene (n = 334,689 gene–subclass combinations) across 26 subclasses, analysis was performed in three steps. (1) Model selection; (2) clustering; and (3) prediction of average trajectory as illustrated in Fig. 2a. The details of these steps are explained below.

Model selection

Our objective was to identify a single optimal model that could capture age-related non-linear trends in gene expression across all 26 subclasses. We achieved this by fitting 12 different models to expression of i subclass and j gene for a total of 334,689 gene–subclass combinations using the dream function as described in equations (5)–(8).

The variable n in equations (6) and (8) represents degrees of freedom, which are 2:6, and the poly function is from the stats (v.4.4.1) R library, which produces orthogonal polynomials with a given degree of freedom in the d.f. argument. To identify the optimal model, we collected the BIC values from each of the 12 models and searched for the minimum values. Although BIC values across all models were not conclusive for glia and other subclasses, a notable decrease in BIC was observed for the poly(log₂(age + 1), d.f. = 2) + covariates model specifically in neurons, particularly INs (Supplementary Fig. 11a). As approximately 60% of the subclasses are neuronal, for simplicity, we adopted the poly(log₂(age + 1), d.f. = 2) + covariates model as the optimal lifespan model for all 26 subclasses. All of the steps in this section were performed using the dream function from the dreamlet (v.1.1.17) R package. The output from the optimal model consists of linear \({\mathrm{coef}}_{{ij}}^{1}\) and non-linear coefficient \({\mathrm{coef}}_{{ij}}^{2}\) from the model for 334,689 gene–subclass combinations.

Clustering

Our goal was to identify the optimal number of unique characteristic curves for the lifespan trends of 334,689 gene–subclass combinations. To achieve this, we applied k-means clustering to \({\mathrm{coef}}_{{ij}}^{1}\) and \({\mathrm{coef}}_{{ij}}^{2}\) from the dreamlet summary statistics of the optimal lifespan model for 334,689 gene–subclass combinations. Clustering all coefficients revealed that k = 10 provided the optimal clusters with non-overlapping lifespan trends (Supplementary Fig. 11b). Supplementary Fig. 11c shows the stratification of all genes across these 10 clusters. For all downstream analyses, we focused on the 135,120 gene–subclass combinations that remained after applying study-wide multiple testing correction (FDR < 0.05) from the optimal model on the 334,689 gene–subclass combinations. Supplementary Fig. 11d illustrates the stratification of lifespan-aDEGs (135,120 genes) across the 10 trajectories. The output table containing \({\mathrm{coef}}_{{ij}}^{1}\) and \({\mathrm{coef}}_{{ij}}^{2}\), along with other summary statistics from the dreamlet tool for each gene and subclass, and trajectory cluster numbers, is provided in Supplementary Data 7.

Prediction of average lifespan trend

Finally, to visualize an average trajectory per cluster, we calculated the mean of the coefficients for all genes within each cluster, as follows:

where N is the total number of genes in the mth cluster and k is the j gene in ith subclass in the mth cluster. Next, we obtained non-linear form of age using \(\hat{\mathrm{Age}} \sim 0+\mathrm{poly}({\log }_{2}({\rm{a}}\mathrm{ge}+1),\text{d.f.}=2)\). Using the mean of the coefficients and the polynomial form of age \(\hat{\mathrm{Age}}\), we obtained lifespan trajectory for m cluster: \(\mathrm{lifespan}\,{\mathrm{trajectory}}_{{m}}\,=\) \(\hat{\mathrm{Age}}\times ({\mathrm{coef}}_{m}^{1},{\mathrm{coef}}_{m}^{2})\), as shown in Fig. 2b. Supplementary Fig. 12b shows predicted average trajectory for each subclass in a trajectory cluster using the mean of subclass coefficients \(\mathrm{lifespan}\,{\mathrm{trajectory}}_{{m},{j}}=\hat{\mathrm{Age}}\times ({\mathrm{coef}}_{m,j}^{1},{\mathrm{coef}}_{m,j}^{2})\). To further obtain biological insights on genes within trajectory clusters, we performed gene set pathway analysis using the enrichr function from the enrichR (v.3.2) R library. The database used was GO Biological Process 2021. Supplementary Fig. 13a shows an example of pathway enrichment from trajectory 1 and 10. The full table of pathways is provided in Supplementary Data 8.

Integrated analysis of lifespan replication dataset

To validate the age-associated transcriptomic changes identified in Fig. 1h, we integrated a human DLPFC snRNA-seq cohort of neurotypical controls from three published datasets^6,7,17, comprising 306 donors spanning gestational week 22 to 101 years of age.

Integration and annotation

Count matrices from multiple snRNA-seq datasets were combined into a single Scanpy object (v.1.9.3) for joint processing and filtering. Genes expressed in fewer than five nuclei across all batches were excluded. Potential doublets were identified and removed using Scrublet (v.0.2.3) with ten principal components per batch. To mitigate sampling bias from variable sequencing depths across studies, we downsampled the integrated data to 1,000 UMI counts per nucleus. Downsampling involved random sampling without replacement, and nuclei with less than 1,000 total UMIs were excluded. The count matrix was then scaled to counts per million after downsampling and transformed using natural log plus one. HVGs (n = 6,000) were selected from protein-coding genes located in the euchromosome using Scanpy’s highly_variable_genes function. Dimensionality reduction was performed by PCA, retaining components explaining 50% of the variance. A neighbourhood graph was constructed using 100 neighbours on the reduced components, followed by generation of a two-dimensional UMAP embedding. Cell type annotation was performed at the class and subclass levels using scvi-tools (v.1.2.0)⁶², transferring labels from our in-house reference cohort.

Replication of aDEGs

To assess the reproducibility of aDEGs identified in the lifespan DLPFC cohort, we applied the dreamlet-enabled limma modelling framework to subclass- and age-group-specific pseudobulk expression profiles in the replication dataset. Our covariate model was expression ~ age + batch + sex, where batch constitutes the dataset resource (Supplementary Fig. 9c). This analysis followed the same modelling strategy as used in the discovery dataset (see the ‘Pseudobulk data aggregation and covariate selection’ section), allowing consistent estimation of age effects across datasets. We next evaluated replication at three levels (Supplementary Fig. 10): (1) using all expressed genes to compute the Spearman correlation of age-associated expression changes between the discovery and replication datasets stratified by subclass and age group; (2) restricting the analysis to aDEGs (FDR < 0.05) identified in the discovery dataset to assess consistency of direction and magnitude of effects; and (3) calculating π1 statistics, which estimate the proportion of discovery aDEGs that show consistent effects in the replication dataset. Notably, for the validation analysis, three fetal donors and their nuclei were excluded to match the age range of the discovery dataset.

Replication of lifespan trajectories

To assess the replication of lifespan transcriptomic trajectories, we implemented two complementary strategies. In the first strategy, we evaluated the replicability of gene-level lifespan trajectories by projecting the age-associated gene expression models derived from the discovery dataset onto the replication dataset. Specifically, for each subclass and gene with a significant age effect (FDR < 0.05) in the discovery cohort, we extracted the age spline coefficients (β₁, β₂) from the dreamlet model and applied them to the transformed ages (log₂(age + 1), d.f. = 2) of the replication donors. We then quantified the agreement between predicted and observed expression levels using R² (squared Pearson correlation) as shown in Supplementary Fig. 15a. As a second, independent strategy, we applied the full trajectory discovery pipeline as described in the ‘Lifespan trends of gene expression’ section to the replication dataset (Supplementary Fig. 15b). This included reaggregating pseudobulk expression data, computing non-linear coefficients using the optimal model and reclustering genes into ten trajectories using the same clustering parameters as in the discovery. We then compared the resulting trajectory shapes and gene memberships between the two datasets to evaluate the concordance of age-associated changes as shown in Supplementary Fig. 15c,d. We also quantified replication performance across trajectories by plotting the proportion of true positives (π₁) recovered in the replication dataset for each trajectory as shown in Supplementary Fig. 15e.

Statistical power analysis

To evaluate the sensitivity of our dataset to detect age-associated changes for each age group and across lifespan, we conducted both theoretical and empirical power assessments.

Theoretical power estimation

We first estimated statistical power using Cohen’s d, sample size and a range of moderated R² values (0.5 to 1), which reflect the proportion of gene expression variance explained by age. Power calculations were Bonferroni corrected for multiple testing, accounting for 26 subclasses and a minimum of around 11,700 expressed genes per subclass across age groups and whole lifespan. As shown in Supplementary Fig. 8a, with a sample size of ≥50, we achieved >80% power to detect effects of d ≥ 1.2 across a range of biologically plausible R² values. Given that all four of our age-defined groups include 50–100 samples—development (n = 53), young adulthood (n = 54), middle adulthood (n = 95) and late adulthood (n = 82)—our design is sufficiently powered to support robust group-level comparisons.

Empirical power evaluation

To assess empirical variation in power across cell subclasses, we considered both technical and biological factors that influence sensitivity to detect age-associated differential expression. First, we examined the relationship between the number of expressed genes and the number of donors per subclass. Supplementary Fig. 8b,e shows that subclasses with greater donor coverage tend to exhibit higher numbers of expressed genes, reflecting improved transcriptome detection with increased sample size, particularly in the late adulthood group. By contrast, the development group displays substantial variability across subclasses, especially among ENs, probably due to underlying biological heterogeneity during early brain development. We further evaluated how these factors impact aDEGs detection sensitivity. As shown in Supplementary Fig. 8c,f, there is a positive association between the number of aDEGs and the number of expressed genes per subclass. Finally, we assessed the relationship between aDEG counts and the magnitude of age-associated effect sizes. Supplementary Fig. 8d,g reveals an inverse relationship: low-powered subclasses (for example, SMC, PC, VLMC) with fewer expressed genes and low aDEG counts show inflated effect sizes, consistent with the winner’s curse, while well-powered subclasses such as neurons and glia exhibit smaller, more stable effect estimates.

Finally, we evaluated power using lifespan aDEGs (defined by F-statistics across both linear and non-linear age terms) across subclasses. Supplementary Fig. 8h displays the overall distribution of donor coverage and gene expression for each subclass using this lifespan non-linear optimal model. Supplementary Fig. 8i,j extends our analysis by correlating lifespan aDEG counts with the number of expressed genes and the mean absolute age effect size.

Degree of sharing across the EN, IN and glial subclasses

To identify genes with similarity in effect sizes across subclasses, the association statistics from dreamlet were not sufficient. For example, if the age-associated effect size for a gene was significant in one subclass but not in another, one might have concluded that the age effect was subclass-specific. However, the absence of a significant effect in another subclass did not necessarily mean that the effect size was zero. This scenario often occurs when statistical power is limited, or varies between subclasses. To overcome this limitation, we used the mashr²⁴ (v.0.2.79) R library, which used an empirical Bayes approach to learn patterns of similarity in effect sizes across subclasses and then leveraged these prior patterns to improve the accuracy of effect-size estimates.

We assessed gene sharing across 9 EN, 7 IN and 4 glial subclasses. The EN group included 9 subclasses (EN_L6_CT, EN_L5_6_NP, EN_L6B, EN_L3_5_IT_1, EN_L3_5_IT_2, EN_L3_5_IT_3, EN_L2_3_IT, EN_L6_IT_1, EN_L6_IT_2), the IN group included 7 subclasses (IN_LAMP5_RELN, IN_LAMP5_LHX6, IN_ADARB2, IN_VIP, IN_PVALB_CHC, IN_PVALB, IN_SST) and the glial group included 4 subclasses (astro, oligo, OPC, micro). Hereafter, we refer to the EN, IN and glial subclasses sets as C_E, C_I and C_G, respectively. In this section, we conducted two analyses: (1) we quantified the degree of sharing for each gene based on the similarity of age-associated effect-size patterns across the EN, IN and glial subclasses using the run_mash function built in dreamlet adapted from mashr; and (2) we calculated the tau score to quantify the cell specificity of each gene, enabling us to stratify shared and non-shared genes.

Mashr analysis

Using the age-associated summary statistics per age group from Supplementary Data 4, we built two matrices with genes as rows and subclasses as columns, containing log₂[FC] (β_j,i) and standard errors (s.e._j,i). The number of genes included was as follows: development (n = 23,914 genes), young adulthood (n = 26,107 genes), middle adulthood (n = 26,809 genes) and late adulthood (n = 24,795 genes). Any gene without summary statistics for a subclass was replaced with 0.

We applied the empirical Bayes framework implemented in mashr to jointly model effect sizes across subclasses, estimating posterior means and variances by combining prior covariance structures with observed sampling variance. This approach improves estimation of true age-associated effects while explicitly modelling sharing patterns across EN, IN and glial subclasses. For each gene–subclass pair, mashr reports the local false sign rate (lfsr), defined as the posterior probability that the true effect has the opposite sign from the estimated effect. We converted lfsr values to posterior probabilities of sign concordance and quantified the degree of sharing within each grouped class (EN, IN, glia) as the product of concordance probabilities across subclasses belonging to that class.

To ensure robustness, we also required that each gene show nominal significance (dreamlet P < 0.05) in at least two subclasses within a grouped class. Composite posterior probabilities were then stratified into ten equally sized bins to visualize sharing patterns, with higher bins indicating stronger cross-subclass concordance.

Using a threshold of composite posterior probability of >0.9, we identified genes with highly shared age-associated effects within each grouped class across four age groups (Supplementary Fig. 16a). The distribution and subclass-level concordance of shared and non-shared genes are shown in Supplementary Fig. 17a,b and Extended Data Fig. 2b. Final composite posterior probabilities for each gene, grouped class and age group are provided in Supplementary Data 11.

Using these shared genes per age group per grouped class, we identified a network of genes encoding significant protein–protein interactions (PPI) using STRING-db⁶³ (v.12.0). Supplementary Fig. 20a–c show the significant PPIs with highest confidence interaction score of >0.9 and genes in networks of at least five genes. The clusters of each PPI per grouped class per age group are given in Supplementary Data 13.

Cell-specificity score

We next reasoned that shared genes with a composite posterior probability of >0.9 have significantly lower cell specificity compared to genes with <0.9. To validate this, we estimated the tau score adapted from GTEx²⁵ studies of each gene for each grouped class per age group. The tau score indicates how specifically a gene is expressed across various subclasses. In other words, genes with a tau score close to 1 are more specifically expressed in one subclass, while those with a tau score closer to 0 are equally expressed across all subclasses within a grouped class. To do this, we integrated pseudobulk gene expression data from a grouped class per age group using the stackedAssays function, resulting in expression of h genes and i subclass × k donors—specifically, 17,205 × 2,393; 17,205 × 1,945 and 17,205 × 1,129 for EN, IN and glia respectively. The analysis was limited to protein-coding genes. Next, we estimated the tau score of each gene across all subclasses. The scores were kept for cell specificity measurement and retained for only those genes with a sum of median expression value across all subclasses of >10 cpm. Supplementary Figs. 18c and 19 show the distribution of tau scores stratified by shared and non-shared genes for each grouped class and each age group for a list of genes from mashr analysis and aDEGs. The tau score of each age group across three grouped classes is provided in Supplementary Table 4. The gene set pathways of shared genes, shown in Extended Data Fig. 2e, were obtained using the enrichr function from the enrichR (v.3.2) R library and the full table is provided in Supplementary Data 12. The analysis used the GO Biological Process 2021 database and was limited to protein-coding genes.

Pseudotime analysis and dynamically expressed gene identification and downstream analysis

For pseudotime analysis, we integrated the lifespan DLPFC dataset (current study) with the published snRNA-seq data spanning gestation to adulthood⁶, according to the data integration procedure described in the ‘Integrated analysis of lifespan replication dataset’ (integration and annotation) section. In contrast to the replication analysis, annotation transfer was not performed; instead, we assessed dataset similarity post-integration. This resulted in a combined dataset of 311 donors. The UMAT embedding was constructed as previously described⁶. Trajectory reconstruction and identification of DEGs along trajectories (traDEGs) was performed as previously described⁴. For each lineage, corresponding cells were selected and genes not observed in ≥5 nuclei were excluded. In total, 6,000 HVGs were selected, data dimensions were reduced using PCA to components explaining 50% of the variance and the UMAT embedding was recalculated. Pseudotime trajectory analysis was then conducted using Monocle3 (v.1.0.0) based on the UMAT embedding. The shortest path between the developmental node and the node in the mature subclass/subtype clusters was isolated as the corresponding trajectory graph. Cells along the trajectory were selected, and traDEGs were identified using Monocle3’s modified graph_test function with Moran’s I test, including covariates (sex, batch, PMI and log_n_genes) to ensure that the results were not affected by uneven contributions from different individuals and conditions. Genes with an adjusted P < 0.05 and Moran’s I ≥ 0.05 were considered statistically significant DEGs. To cluster the DEGs in each lineage, single-cell expression data along each trajectory were compressed using a sliding window along pseudotime, averaging the expression of neighbouring cells to generate 500 metacells per trajectory. Each gene’s expression was then modelled using a generalized linear model (\(\mathrm{expression} \sim \mathrm{splines}::\mathrm{ns}(\mathrm{pseudotime},\text{d.f.}=3)\)), and k-means clustering was performed on the fitted expressions.

GO-term analysis of traDEG clusters was performed using enrichR (v.3.2) and the GO Biological Process 2023 dataset. To investigate whether the traDEG clusters have a role in neurological, psychiatric and other traits, we quantified their colocalization with common risk variants from 24 GWAS (Supplementary Table 2) using MAGMA analysis. Single-cell disease scores along each trajectory were compressed into 500 metacells per trajectory. Those scores were then modelled using a generalized linear model \(\mathrm{diseasescore} \sim \mathrm{splines}::\mathrm{ns}\) (pseudotime, d.f. = 3). In contrast to MAGMA, which assesses gene-level GWAS enrichment across predefined cell types, scDRS provides per-nucleus scores by weighting gene expression with GWAS-derived effect sizes.

Validation of traDEG modules with spatial transcriptomics

To validate spatial expression patterns of traDEG modules, we analysed Visium spatial transcriptomics data from 10 adult donors (aged 33–62 years; three cortical sections per donor, 30 sections total) from a publicly available dataset³³. Data were accessed through the spatialLIBD R package (v.1.19.11)⁶⁴ using fetch_data(type = “spatialDLPFC_Visium”).

For each traDEG module, multi-gene expression was summarized per Visium spot using PCA implemented in vis_gene with multi_gene_method =”pca”. The first principal component (PC1) was used as a composite score representing module expression. Region-specific enrichment was visualized by computing average PC1 scores per BayesSpace-defined spatial cluster (k = 9).

Spatial validation was performed for all glial and neuronal modules. Visualization outputs included heat maps of average PC1 scores across donors, PC1 score maps for representative sections (Br2743_mid) and violin plots depicting variability across all sections.

Enrichment of brain- and non-brain-related risk genes in age and pseudotime associated genes

All analyses to evaluate the enrichment of brain and non-brain related traits were conducted using MAGMA (v.1.08b)⁶¹. MAGMA provided robust gene-set enrichment statistics, reinforcing the biological relevance of our findings in the context of ageing and disease; it calculates gene-level P values for each gene and trait by assessing the joint association of all single-nucleotide polymorphisms within the gene region, while accounting for linkage disequilibrium (LD) between single-nucleotide polymorphisms. Gene regions were defined with a window of 35 kb upstream and 10 kb downstream, and linkage disequilibrium estimates were derived from the European panel of the 1000 Genomes Project⁶⁵ (phase 3). MAGMA then applies a linear regression framework to determine whether DEGs show stronger associations with GWAS traits compared with the rest of the genome. Genes overlapping the MHC region (chromosome 6: 25–35 Mb) were excluded from the analysis. Heat maps of brain- and non-brain-related traits were generated using the MAGMA pipeline as described. Outputs are provided in the Supplementary Data 6 (Fig. 1j and Supplementary Fig. 7g), Supplementary Data 9 (Fig. 2e,g and Supplementary Fig. 12c) and Supplementary Data 15 (Figs. 3f and 4c,g and Extended Data Fig. 5e,i). For Supplementary Data 15, we extracted significant genes (FDR < 0.05) from MAGMA gene-level outputs for each trait and overlapped them with pseudotime-inferred gene modules to enable detailed gene-level exploration of trait enrichment.

Inference of transcriptional regulators of pseudotime gene modules

To identify TFs regulating pseudotime-driven gene modules, we used the ULM approach in decoupleR (v.2.9.7)³¹. This method fits a linear model for each TF, using TF–target interaction weights from the regulatory network CollecTRI⁶⁶, and calculates a t-statistic as a measure of TF enrichment. Only TFs with a minimum of 3 target genes for modules containing fewer than 50 genes, or 5 target genes for larger modules, were considered. TF activity inference was performed independently for each lineage-specific pseudotime trajectory. TFs were considered significant if they showed enrichment (P < 0.05) in at least 2% of pseudotime bins along a trajectory, and significant TFs were aggregated across lineages for downstream analyses.

Protein–protein interaction networks were constructed using STRING-db (as described in the ‘Degree of sharing across the EN, IN and glial subclasses’ section) by combining genes from each traDEG module with their corresponding decoupleR-inferred TFs. To prioritize key regulators, module-specific networks were filtered to retain TFs within connected components of size ≥ 5 and node degree ≥ 3. The resulting set of 44 key TFs was further evaluated for biological relevance by comparison to 582 enhancer-driven regulons (eRegulons) identified with SCENIC+ in the developing human neocortex³².

Rhythmicity analysis

Before rhythmicity analysis, the TOD for each individual was normalized to a ZT scale as described previously⁶⁷. In brief, the TOD for each individual was collected at local time then converted to coordinated universal time by adjusting for time zone and daylight saving time. Coordinated universal time was further adjusted to account for the longitude and latitude of death place. Each individual’s TOD was then set as ZT = t h after previous (if t ≥ 18) or before next (if t ≥ −6) sunrise. The distribution of participant TODs across age groups is shown in Extended Data Fig. 6a. The infancy, childhood and adolescent groups had too few participants with known TODs to run further rhythmicity analyses, and the young adulthood group had large gaps in their distribution. As such, we excluded the infancy, childhood and adolescent groups and combined the young adulthood + middle adulthood groups (n = 119), which was then compared with the late adulthood group (LA; n = 77) (Fig. 6a).

To measure the rhythmicity in i subclass and j genes and for other downstream analyses in this section, we first regressed out the effects of identified technical and biological covariates in the previous section along with age using equation (9) as described below.

Then, using the DiffCircaPipeline (v.0.0.0.9000) workflow⁶⁸, 24 h rhythms in gene expression were detected and compared between the young adulthood + middle adulthood and late adulthood groups. First, samples were ordered by TOD and residualized expression for each transcript was fit to a cosinor model separately in each group. The most robust rhythmic genes were identified as those that survived multiple comparison testing (FDR-corrected P < 0.05, Benjamini–Hochberg method). For quality control, we also performed rhythmicity analyses on 100 permutations of our data with scrambled TOD values within both the young adulthood + middle adulthood and late adulthood groups (Supplementary Data 18). We then used these permutations to calculate empirical P values: [(number of permutation R² values ≥ observed R² value) + 1]/[(number of permutations) + 1].

For a deeper analysis, we used a P < 0.01 threshold within DiffCircaPipeline to categorize the type of rhythmicity (rhythmic in young adulthood + middle adulthood, rhythmic in late adulthood, both or arrhythmic) of each gene (Supplementary Data 17). These categories were then used to determine which genes to perform differential goodness of fit and parameter tests on to reduce the likelihood of a type I error⁶⁸. Specifically, differences in goodness of fit (ΔR²) between the young adulthood + middle adulthood and late adulthood groups were determined through a permutation test (1,000 permutations) in transcripts identified as rhythmic in young adulthood + middle adulthood, late adulthood or both (Supplementary Data 20). Moreover, a global differential parameter test was performed in transcripts identified as rhythmic in both groups, followed by post hoc analyses to determine whether differences were due to the rhythmic parameters amplitude, MESOR or phase (Supplementary Data 19). enrichR was then used to identify pathways enriched in rhythmic genes (Supplementary Data 21). This was done for all rhythmic genes within a subclass in young adulthood + middle adulthood and late adulthood and then separately for rhythmic genes that peaked between 0–8 h ZT and 20–24 h ZT (the sunrise group) and 8–20 h ZT (the sunset group).

Ethics statement

All procedures and research protocols were approved by the Institutional Review Boards (IRBs) of the collaborating institutions: the Ethics Committee at Icahn School of Medicine at Mount Sinai, James J. Peters Department of Veterans Affairs Medical Center and the Human Brain Collection Core at the National Institute of Mental Health (NIMH). All research was performed in accordance with relevant guidelines and regulations. Informed consent was obtained from all participants.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.