Lifespan single-cell transcriptomic atlas of the human prefrontal cortex

The organization of molecular and cellular programs across the human lifespan remains a central challenge in understanding brain development, ageing and disease susceptibility. Although single-cell transcriptomic studies^4,5,6,7 have transformed our understanding of cellular diversity in the human brain, how transcriptional programs evolve continuously from early development through late adulthood within defined cortical regions remains poorly characterized. The dorsolateral prefrontal cortex (DLPFC), which supports higher cognitive functions³, is particularly well suited for investigating lifespan dynamics given its prolonged maturation and sensitivity to age-related decline.

To address this gap, we generated a single-nucleus transcriptomic atlas of DLPFC spanning the full lifespan, profiling over 1.3 million nuclei from neurotypical donors. This resource, produced as part of the PsychAD Consortium^8,9, enables systematic analysis of cell-type-specific transcriptional programs across development, adulthood and ageing. Using integrated differential expression, trajectory modelling and spatial transcriptomic analyses, we identify three major phases of molecular change, define age-dependent gene modules within neuronal and glial lineages, and link these programs to genetic risk for psychiatric and neurodegenerative disorders. We further identify age-related reorganization of circadian gene expression, marked by loss of neuronal clock synchronization and the emergence of glial stress-associated rhythms. Together, these findings establish a framework for understanding how cellular programs transition from resilience to vulnerability across the human lifespan.

To construct a lifespan-resolved transcriptomic atlas of the DLPFC at the single-nucleus resolution, we profiled single-nucleus RNA-sequencing (snRNA-seq) libraries from 284 neurotypical donors spanning infancy (0–1 years old), childhood (2–11 years old), adolescence (12–19 years old), and young (20–39 years old), middle (40–59 years old) and late (≥60 years old) adulthood^8,9 (Fig. 1a, Supplementary Fig. 1a, Supplementary Table 1, Supplementary Data 1 and Supplementary Notes). For downstream analyses (Fig. 1b), infancy, childhood and adolescence were consolidated into a single development group to increase the nucleus counts for the early stage, therefore yielding four broad age groups. After stringent quality-control filtering (Supplementary Fig. 2a and Methods), we retained 1,307,674 single-nucleus transcriptomes. These were annotated into eight major cell classes: excitatory (EN) and inhibitory (IN) neurons, astrocytes (astro), immune cells (immune), mural cells (mural), endothelial cells (endo), oligodendrocytes (oligo) and oligodendrocyte progenitor cells (OPC). These classes were further resolved into 27 subclasses and 65 subtypes (Fig. 1a and Supplementary Fig. 1b). Subclass identities were validated using canonical marker genes and alignment with external references, including the dataset from a previous study¹⁰ (Supplementary Fig. 3).

To quantify age-related transcriptional variation across subclasses and age groups, we applied a linear mixed model using the variancePartition¹¹ tool (Supplementary Fig. 2b,c) on pseudobulk gene expression aggregated at the donor-by-subclass level. Age explained the greatest transcriptional variance during development across all subclasses (mean, 4%), reflecting dynamic transcriptomic shifts associated with postnatal neurogenesis and gliogenesis (Fig. 1c and Supplementary Fig. 2d). The second-largest age-associated effect was observed during late adulthood (mean, 2%) (Supplementary Fig. 2d and Supplementary Data 3), consistent with emerging transcriptional alterations linked to ageing¹².

We next investigated the enrichment of subclass-level transcriptomes for traits relevant to brain health (Methods and Supplementary Fig. 4). Among 24 traits grouped into psychiatric, neurological and other (metabolic and immunological) categories (Supplementary Table 2 and Supplementary Data 2), psychiatric traits were predominantly enriched in neuronal subclasses. Major depressive disorder, autism spectrum disorder, attention deficit hyperactivity disorder and educational attainment also demonstrated significant enrichment in OPCs. By contrast, neurological and other traits were most strongly associated with immune, mural and astrocytic subclasses, with immune cells (adaptive, microglia (micro) and perivascular macrophages (PVM)) demonstrating a particularly pronounced effect¹³. Notably, obesity showed enrichment in both the EN and IN classes, highlighting a shared genetic architecture between metabolic and psychiatric traits¹⁴. Together, these analyses reveal how transcriptional programs in distinct DLPFC cell types relate to genetic risk for psychiatric, neurological and metabolic disorders, establishing this lifespan atlas as a framework for investigating disease vulnerability and cellular resilience.

Nucleus composition analysis identified prominent shifts in the EN and IN subclasses during development (infancy–adolescence; Fig. 1d), with peak age-related variance at 9.7% versus 0.22–1.50% in adulthood (Fig. 1e and Methods). Modelling abundance changes using crumblr¹⁵ (adjusted for sex, post-mortem interval and sample source; Supplementary Fig. 5a) identified development and young adulthood as periods of maximal compositional change (Fig. 1f and Supplementary Table 3b–e). Notably, a transition point emerged at 24 years of age, after which compositional changes declined sharply, marking a shift from active developmental remodelling to relative cellular stability in the DLPFC (Methods and Supplementary Fig. 5c,d). While most subclasses stabilized in adulthood, the SST-positive IN (IN_SST) and OPC subclasses continued to exhibit significant abundance shifts in late adulthood, indicating that specific neuronal and glial populations retain compositional plasticity throughout ageing.

Lifespan trajectories of nucleus composition fit logarithmic rather than linear age models for nearly all subclasses (Supplementary Figs. 5b and 6), revealing two major patterns: log-increasing trajectories (ENs, astro, oligo) and log-decreasing trajectories (IN, OPC, micro; Extended Data Fig. 1a). These non-linear trends mirror key principles of developmental biology, in which rapid postnatal proliferation and synaptic remodelling drive accelerated compositional change, followed by stabilization after early adulthood (Extended Data Fig. 1b). Across 26 subclasses, 24 exhibited significant compositional shifts after accounting for cellular hierarchy (Extended Data Fig. 1c and Supplementary Table 3a). To validate these findings, we performed RNAscope assays targeting SST (IN_SST) and GPR37 (oligo) across representative age groups spanning development (4–9 months, 9 years), young adulthood (38 years), middle adulthood (51–53 years) and late adulthood (65–86 years). These assays revealed a progressive decline in SST⁺ interneurons and a marked increase in GPR37⁺ oligodendrocytes with age (Fig. 1g and Extended Data Fig. 1d,e), consistent with predictions from the logarithmic models.

Pathway enrichment analysis revealed distinct molecular programs along compositional trajectories (Extended Data Fig. 1f). Log-increasing subclasses were enriched for neuronal development and synaptic organization pathways, with ENs showing particularly strong enrichment for synaptic organization and neurogenesis. INs showed lower developmental pathway enrichment but substantial enrichment for neurotransmitter metabolic and biosynthetic processes, consistent with specialized inhibitory circuit roles. Log-decreasing glial subclasses were enriched for proliferation, apoptotic and neuroinflammatory processes. Together, coordinated molecular programs of cellular expansion and attrition orchestrate DLPFC remodelling, with developmental windows characterized by intense remodelling and age-related stabilization.

The age-related variance in subclass transcriptomes (Fig. 1c) suggested dynamic transcriptional remodelling across the lifespan. Using dreamlet¹⁶ at the pseudobulk level (Methods), we identified three distinct transcriptional phases in the DLPFC (Supplementary Fig. 7a,b and Supplementary Data 4). The first phase, spanning development, exhibited widespread remodelling with 8,223 age-associated differentially expressed genes (aDEGs) dominated by neuronal subclasses (41.2% EN; 40.5% IN) (Fig. 1h). The second phase, encompassing young and middle adulthood, showed transcriptomic stability with only 27 and 1 aDEGs, respectively. By contrast, the third phase during late adulthood showed renewed remodelling, primarily affecting glia (426 out of 735 aDEGs) versus neurons (246 out of 735) (Fig. 1h). Pronounced age-associated effect sizes during development compared with during late adulthood further highlighted magnitude of early-life changes (Supplementary Fig. 7c). Power analyses indicated that the study design (at least 50 donors per age group; around 11,700 median expressed genes) was sufficient to detect age-associated effects across all groups (Supplementary Fig. 8a–g), supporting the robustness of the observed triphasic pattern. Consistent with this, independent validation using a lifespan replication dataset comprising 306 neurotypical donors (gestational week 22 to 101 years; 2,039,078 nuclei) from published sources^6,7,17 confirmed reproducible aDEGs with concordant effect directions and magnitudes, particularly during development and late adulthood (Methods and Supplementary Figs. 9 and 10).

Developmental aDEGs showed distinct functional signatures. Downregulated aDEGs were enriched for neurodevelopmental processes (axon guidance, neuronal differentiation; Fig. 1i), whereas upregulated aDEGs were associated with immune and metabolic pathways across the EN_L6B, oligo, adaptive, micro and endo subclasses (Fig. 1i). Psychiatric risk genes for schizophrenia and major depressive disorder were enriched among downregulated aDEGs but not upregulated aDEGs¹⁸ (Fig. 1j and Supplementary Data 6). Downregulated aDEGs showed higher loss-of-function intolerance (probability of loss-of-function intolerance (pLI) scores¹⁹) compared with upregulated aDEGs (Supplementary Fig. 7d), indicating stronger purifying selection. In late adulthood, glial upregulated aDEGs were enriched for leucine metabolism, DNA binding and protein localization, whereas downregulated aDEGs in the OPC and IN subclasses were associated with synaptic transmission (Supplementary Fig. 7f and Supplementary Data 5). Glial and IN_ADARB2 aDEGs showed significant enrichment for Alzheimer’s disease (AD) and other neurological traits (Supplementary Fig. 7g and Supplementary Data 6), with late-adulthood aDEGs showing no significant pLI association except in OPCs (Supplementary Fig. 7e).

Together, these analyses reveal temporally segregated transcriptional programs in the human DLPFC characterized by neuronal specification under strong evolutionary constraint during development, relative transcriptional stability in midlife, and late-life reactivation of glial networks linked to immune regulation and neurodegeneration. This triphasic framework establishes a molecular foundation for understanding age-related transcriptional shifts in cortical maturation, cognitive function and disease vulnerability.

After observing significant transcriptomic changes during development and late adulthood, we next examined temporal trajectories of aDEGs to characterize their expression dynamics across the lifespan. Following the workflow illustrated in Fig. 2a, we evaluated 12 alternative ageing models (linear, log-transformed and polynomial forms) using subclass-level pseudobulk matrices (334,689 gene–subclass combinations; median, 13,227 genes per subclass). Based on the median Bayesian information criterion (BIC) across subclasses, we selected the second-degree polynomial of log-transformed age as optimal, providing the best fit, particularly in neuronal subclasses (Methods and Supplementary Fig. 11a). Power analysis confirmed that this model was sufficiently powered to detect non-linear age-associated effects (Supplementary Fig. 8h–j).

Approximately 40.3% (135,120 out of 334,689) of gene–subclass combinations showed significant age-related changes, which clustered into ten distinct trajectories (Fig. 2b, Supplementary Fig. 11b–d, Supplementary Fig. 12a,b and Supplementary Data 7). Trajectory 1 (16,188 genes from 23 subclasses) captured basic cellular processes, while trajectory 5 (5,172 genes from 26 subclasses) was enriched for WNT signalling, extracellular matrix organization and synaptic transmission across glial and neuronal subclasses (Supplementary Fig. 13a and Supplementary Data 8). These findings reveal divergent molecular programs associated with age-related dynamics across the lifespan.

Developmental aDEGs mapped predominantly to trajectories 2 and 7 in the EN and IN subclasses (Fig. 2c). Trajectory 2, containing downregulated developmental aDEGs, was enriched for synaptic signalling, glutamatergic transmission and CNS development (Fig. 2d). Notably, genes within this trajectory showed significant overlap with brain-related genome-wide association study (GWAS) traits including schizophrenia, bipolar disorder and educational attainment (Fig. 2e, Supplementary Fig. 12c and Supplementary Data 9). Most developmental aDEGs (67.6% across trajectories 1, 2, 6 and 7) peaked at around adolescence (around 12–13 years old) and remained stable, suggesting that early transcriptomic remodelling establishes neuronal foundations for later brain function and psychiatric vulnerability.

Late-adulthood aDEGs concentrated in trajectories 3–5 and 8–10 (Fig. 2f) exhibited three patterns: pronounced development-to-midlife stabilization with late-life reactivation (trajectories 3 and 8), progressive lifespan change (trajectories 5 and 10) and reversal after early peaks (trajectories 4 and 9) (Supplementary Fig. 13b). Trajectories with upward trends (4, 8 and 10) were enriched for neurological and immune-related disorder risk genes (Fig. 2g). Trajectory 10 specifically highlighted glial immune activation and interferon signalling in late-life neuroinflammation (Fig. 2h). Microglial genes in trajectory 10 showed enrichment for protein localization and calcium homeostasis, with significant AD risk gene overlap (Supplementary Fig. 14a–c). Moreover, microglial trajectories 3 (chromatin modification) and 6 (autophagy), characterized by age-related decline, also showed AD risk enrichment, suggesting that these programs contribute to AD vulnerability independent of late-life differential expression (Supplementary Fig. 14d and Supplementary Data 10). Among neuronal aDEGs, only IN_ADARB2 showed significant late-adulthood pathway enrichment, with downregulation during both development and late adulthood, indicating diminished synaptic plasticity, a hallmark of neuronal ageing^20,21 (Supplementary Fig. 14e,f and Supplementary Data 10).

To assess trajectory reproducibility, we applied non-linear age coefficients to predict pseudobulk expression in the independent lifespan replication dataset (Fig. 2i and Supplementary Fig. 15). Cross-dataset prediction yielded high concordance, with trajectories 5 and 10 exhibiting the highest R² values in both datasets, indicating robust generalizability. Notably, reproducible trajectories were highly cell-type specific and depleted for housekeeping genes, whereas trajectory 1 (the least reproducible) was enriched for housekeeping genes and least cell-type specific (Fig. 2j). These findings suggest that cell-type-specific lifespan programs are consistently recapitulated across individuals, while housekeeping-dominated trajectories capture greater interindividual variability, highlighting the robustness of specific non-linear trajectories and their cell-type specificity in defining individual-level DLPFC variability.

Age-related transcriptional convergence has been described in bulk human brain tissue²² and in mouse models²³, but its cell-type-specific dynamics in the human DLPFC remain poorly characterized. To address this, we applied mashr²⁴—an empirical Bayes framework that models shared effect sizes across conditions. For each age group, we estimated composite probabilities (P_E, P_I, P_G) summarizing shared age-associated changes across 9 EN, 7 IN and 4 glial subclasses (Extended Data Fig. 2a). Analyses focused on genes with nominally significant age effects in at least two subclasses (Supplementary Fig. 16a and Supplementary Data 11).

The EN subclass showed a marked increase in gene sharing with age, with P_E > 0.9, rising substantially from development to late adulthood (Extended Data Fig. 2b), while the IN and glial subclasses exhibited more modest trends (Supplementary Fig. 17a). Low-sharing genes (P_E < 0.01), such as KCNH4, exhibited divergent age effects (Supplementary Fig. 17b), whereas highly shared genes such as PARP2 (P_E = 0.99) showed consistent changes across all EN subclasses (Extended Data Fig. 2c). Averaging across age groups revealed progressive EN convergence from development through late adulthood (Extended Data Fig. 2d), with moderate changes in the IN and glial subclasses (Supplementary Figs. 17c and 18a,b).

Highly shared genes (P_E, P_I, P_G > 0.9) exhibited lower tau scores²⁵ compared with low-sharing genes, indicating reduced cell-type specificity (Supplementary Fig. 18c, Supplementary Fig. 19 and Supplementary Table 4), suggesting that convergence preferentially involves genes with ubiquitous cellular maintenance roles. Functionally, shared genes in early life were enriched for neuronal adhesion and neurotransmitter transport (Supplementary Data 12), whereas late-adulthood shared genes were enriched for DNA repair (EN), ubiquitination (IN) and antigen processing (glial) (Extended Data Fig. 2e).

Protein–protein interaction networks (Methods) revealed a transition from divergent programs to shared ageing pathways. Young-adulthood EN networks highlighted mitochondrial and axonal genes; IN networks featured synaptic proteins (Supplementary Data 13). By late adulthood, EN and IN networks showed clusters in DNA repair and RNA processing, while glial networks exhibited apoptosis and MHC class I signalling (Extended Data Fig. 2f-h and Supplementary Fig. 20). Most genes in these networks were upregulated with age, consistent with increased demands for genomic stability and immune surveillance. Together, these findings reveal progressive transcriptional convergence in the human DLPFC with age, most pronounced in ENs.

Recognizing the variance observed during development and late adulthood (Fig. 1c and Extended Data Fig. 2b), we performed pseudotime trajectory analysis to delineate cellular dynamics of DLPFC lineages across the lifespan. To capture early developmental transitions, we integrated our dataset with published snRNA-seq spanning gestation to adulthood⁶, resulting in a combined dataset of 1,454,617 nuclei from 311 individuals (Methods and Supplementary Fig. 21a–d). We used uniform manifold approximation and projection (UMAP) of maturation (UMAT)⁶, which restricts neighbour selection to adjacent stages, thereby enhancing the precision of representing transitional processes across the lifespan.

In astrocytes, we identified two distinct age-related patterns (Fig. 3a) corresponding to fibrous astrocytes (FAs; high GFAP expression) and protoplasmic astrocytes (PAs; low GFAP, high SLC1A2) (Extended Data Fig. 3a). Pseudotime trajectories revealed that PAs matured later than FAs (Fig. 3b,c), reflecting their distinct microenvironments: FAs in the white matter interact with myelinated axons and oligodendrocytes, whereas PAs in the grey matter engage with neurons and synapses²⁶. Notably, PAs exhibited consistently higher single-cell disease relevance scores (scDRS) for migraines across pseudotime, especially during maturation and ageing (Fig. 3d). These findings underscore the importance of cellular interactions and temporal dynamics in understanding astrocyte functions and their role in disease susceptibility.

To resolve molecular programs shaping astrocyte trajectories, we identified 827 DEGs along pseudotime trajectories (traDEGs; false-discovery rate (FDR)-adjusted P < 0.05, Moran’s I ≥ 0.05) clustered into six modules reflecting developmental (dev), mature (mat) and ageing processes (Fig. 3e, Methods, Extended Data Fig. 3b and Supplementary Data 14). The dev module was enriched for psychiatric traits including bipolar disorder and educational attainment (Fig. 3f, Methods and Supplementary Data 15), and functional analysis revealed involvement in nervous system development and axonogenesis (Fig. 3g and Supplementary Data 16), emphasizing astrocytic contributions to cortical development and mental health²⁷. By contrast, the ageing module was enriched for neurological and immunological traits and linked to macromolecule biosynthesis and cytokine response, implicating astrocytic programs in neuroinflammation and late-life vulnerability²⁸ (Fig. 3f,g). The PA-specific module showed strong enrichment for schizophrenia and migraine risk genes (Fig. 3f), consistent with the scDRS results (Fig. 3d), suggesting that functional interactions between neurons and PA may contribute to disease aetiology^17,29. Moreover, the FA-specific mat-ageing module was linked to memory, negative regulation of amyloid fibril formation (Fig. 3g) and cellular response to organic cyclic compounds, underscoring their involvement in cognitive support, neuroprotection and cellular defence mechanisms³⁰.

To elucidate regulatory mechanisms shaping these astrocytic programs, we inferred transcription factor (TF) activity for each traDEG module using the univariate linear model (ULM) in decoupleR³¹ (Methods). Candidate TFs were integrated with module genes into PPI networks constructed with STRING-db, and prioritized based on their presence in connected components (size ≥ 5) and high node degree (≥3). Among astrocytic modules, robust TF–gene networks were identified for the dev, FA-specific mat-ageing and ageing modules (Extended Data Fig. 3c), reflecting stage- and lineage-specific regulatory architectures. For example, AR emerged as a key developmental regulator, coordinating a network linked to protein folding (Fig. 3g and Extended Data Fig. 3c). Validation against enhancer-driven regulons (eRegulon) from SCENIC+ in the developing neocortex³² further supported AR’s role, along with additional regulators identified in astrocytic modules (Extended Data Fig. 3d,e). These findings reveal dynamic, stage-specific TF networks that coordinate astrocytic contributions to brain development and disease vulnerability.

To validate the anatomical relevance of pseudotime modules, we analysed Visium spatial transcriptomics data from the human DLPFC³³ (Methods). This approach enabled us to assess whether distinct pseudotime modules exhibit laminar or region-specific enrichment, thereby linking temporal transcriptional trajectories to spatially organized astrocyte functions in situ. PA-specific genes were broadly distributed across cortical layers (Fig. 3h and Extended Data Fig. 3f), aligning with neuronal support functions and psychiatric trait associations (Fig. 3f,g). By contrast, FA-specific genes were enriched in white matter (Fig. 3h and Extended Data Fig. 3f), consistent with neuroprotective roles (Fig. 3g). Together, these findings reveal temporally and spatially resolved astrocytic programs, underscoring subtype- and stage-specific roles in shaping brain health and disease vulnerability across the human lifespan.

We next examined the lifespan dynamics of other glial lineages, focusing on microglia and oligodendrocytes. For microglia, pseudotime trajectory analysis revealed a single developmental path (Fig. 4a and Extended Data Fig. 4a). Trait enrichment analysis showed significant and stable enrichment in several neurological traits, including AD and multiple sclerosis, as well as immunological traits such as rheumatoid arthritis, inflammatory bowel disease and ulcerative colitis (Extended Data Fig. 4b). Similar to astrocytes, developmental genes in microglia were associated with psychiatric traits and neuronal development processes, while ageing-associated genes were linked to immunological functions (Fig. 4b–d, Extended Data Fig. 4c and Supplementary Data 14–16). These findings highlight the dual roles of microglia in supporting brain development and contributing to neurodegeneration and immune-related conditions across the lifespan^34,35. Notably, the progression from developmental to ageing modules in microglia revealed genetic associations with AD that emerge later in life, offering insights beyond those identified in previous GWAS focused on aged donors³⁶ (Fig. 4b–d).

For the oligodendrocyte lineage, we identified a trajectory from OPCs to mature oligodendrocytes (Fig. 4e–h and Extended Data Fig. 4a). Trait enrichment along this trajectory revealed distinct transitions: OPCs showed significant associations with psychiatric traits and obesity, while mature oligodendrocytes exhibited enrichment for AD (Fig. 4g). Moreover, developmental genes in oligodendrocytes were associated with multiple psychiatric disorders and synaptic transmission, whereas ageing genes were linked to myelination, protein localization to plasma membrane and obsessive compulsive disorder (Fig. 4f–h and Supplementary Data 14–16). These patterns highlight the stage-specific functions of oligodendrocytes in maintaining neural circuitry and their roles in lifespan vulnerability to psychiatric and neurological traits.

TF–gene network analysis identified CTNNB1, PPARGC1A and ESR1 as key regulators within the oligodendrocyte developmental module, orchestrating networks tied to synaptic signalling and axon development (Fig. 4h and Supplementary Fig. 22). In microglia, the ageing module featured HDAC4-, RUNX1- and HIF1A-centred networks associated with immune response pathways (Fig. 4d and Extended Data Fig. 4d). Spatial transcriptomic validation revealed that developmental modules in both lineages were broadly distributed across cortical layers, while ageing modules were enriched in white matter regions (Fig. 4i,j and Extended Data Fig. 4e). Together, these findings underscore the temporally and spatially distinct regulatory programs of microglia and oligodendrocytes, linking their lineage-specific dynamics to lifespan disease susceptibility.

We next focused on neuronal lineages, including ENs and INs (Fig. 5 and Extended Data Fig. 5). For ENs, pseudotime trajectory analysis revealed nine developmental paths grouped into three categories based on cortical distribution: upper-layer intratelencephalic projection neurons (upper-IT), deep-layer intratelencephalic projection neurons (deep-IT) and deep-layer non-intratelencephalic projection neurons (deep-non-IT) (Fig. 5a and Extended Data Fig. 5a). These categories displayed distinct temporal patterns, with deep-non-IT ENs maturing earliest, followed by deep-IT and upper-IT ENs (Fig. 5a,b and Extended Data Fig. 5b), consistent with the canonical inside-out pattern of cortical development³⁷.

Across these trajectories, we identified 2,686 genes exhibiting dynamic expression, which clustered into eight distinct modules (Fig. 5c, Extended Data Fig. 5c and Supplementary Data 14). Spatial transcriptomics revealed marked laminar enrichment of these modules: upper-layer-specific modules were localized predominantly to cortical layers II–IV, while deep-layer-specific modules were enriched in layers V–VI (Fig. 5d and Extended Data Fig. 5d). This laminar segregation underscores the spatial distinction of transcriptional programs across EN subclasses and highlights their developmental divergence and functional specializations within cortical circuits.

Functional enrichment analysis revealed that developmental genes were associated with neuron generation and development, while mature and ageing genes were enriched for synaptic transmission and ion transport (Fig. 5e and Supplementary Data 15), essential for synaptic function and structural integrity throughout life³⁸. Deep-layer modules included developmental and maturation programs enriched in neuron projection development and protein folding, with cholesterol metabolism uniquely regulated by SREBF2-centred networks (Fig. 5e and Supplementary Fig. 23). Ageing modules encompassed processes such as extracellular matrix organization and circadian entrainment, with key regulators including NR3C1. Across upper- and deep-layer modules, shared regulators such as CTNNB1, MYC and ESR1 orchestrated these programs, reflecting both common and trajectory-specific transcriptional control. Disease enrichment analysis revealed robust associations between ENs and psychiatric traits (Extended Data Fig. 5e and Supplementary Data 16). Schizophrenia, for example, showed significant associations during developmental and mature stages, suggesting critical windows of vulnerability³⁹. Tourette’s syndrome and obesity exhibited progressive associations during EN maturation, indicating roles for neuronal development in these traits. These results extend previous findings that primarily emphasized IN involvement^40,41. Together, these findings underscore how trajectory-dependent transcriptional programs shape EN development, synaptic maintenance and disease associations.

For INs, pseudotime analysis identified 11 trajectories, grouped into 2 categories based on origin: medial ganglionic eminence (MGE)-drived and caudal ganglionic eminence (CGE)-derived INs (Fig. 5f and Extended Data Fig. 5f). While overall maturation patterns of CGE- and MGE-derived INs were similar, CGE-derived INs exhibited greater variability across trajectories (Fig. 5g and Extended Data Fig. 5g), suggesting distinct temporal dynamics across these IN categories. We identified 2,036 traDEGs, forming eight distinct modules (Fig. 5h, Extended Data Fig. 5h and Supplementary Data 14). Developmental modules were enriched for genes associated with nervous system development, axonogenesis and axon guidance (Fig. 5i and Supplementary Data 15). Moreover, CGE-specific developmental genes were involved in regulation of intracellular signal transduction. Mature and ageing genes were enriched in synaptic transmission and ion transport, contributing to the stability and adaptability of neural circuits throughout life³⁸. Specifically, sporadic trajectories in the mature and ageing modules showed enrichment for extracellular matrix maintenance, echoing patterns observed in L5_6_NP ENs (Fig. 5i,e). TF–gene network analysis revealed modular architectures with distinct functional enrichments, reflecting stage- and subtype-specific regulation (Fig. 5i and Supplementary Fig. 24). Furthermore, enrichment in psychiatric disorders was also observed for IN traDEGs (Extended Data Fig. 5i and Supplementary Data 16), with the enrichment being more pronounced in developmental genes. Notably, PVALB_CHC-specific mature and ageing genes showed significant enrichment in alcoholism and stroke, potentially due to their involvement in response to calcium ions (Fig. 5i). Together, these findings reveal distinct temporal dynamics and disease associations within IN lineages, underscoring their contributions to both developmental and late-life brain disorders.

Given the enrichment of TF networks for circadian entrainment pathways and the established link between circadian rhythm disruption and ageing^42,43, we investigated how age influences circadian transcriptional programs in the DLPFC. Using a cosinor model, we identified 24 h gene expression rhythms in individuals with known times of death (TOD) (Fig. 6a, Methods and Extended Data Fig. 6). This enabled us to compare rhythmicity in young and middle adulthood (n = 119) to late adulthood (LA, n = 77) (Fig. 6a).

In the young adulthood + middle adulthood group, we detected 64 significant rhythmic genes (FDR < 0.05) across 11 subclasses, with 69% (44 out of 64) found in upper-layer IT ENs (Fig. 6b and Supplementary Data 17). Many of these genes are associated with the circadian molecular clock (Extended Data Fig. 6b). By contrast, only VCAN exhibited significant rhythmicity in late adulthood (Fig. 6b). Permutation testing confirmed the robustness of these rhythmic signatures (empirical P = 0.0099; Supplementary Fig. 25 and Supplementary Data 18). Applying a relaxed threshold (P < 0.01) revealed more rhythmic genes in young adulthood + middle adulthood than late adulthood across most subclasses, except for EN_L6_CT, EN_L6B, IN_LAMP5_LHX6, micro and oligo, in which rhythmic genes were more abundant in late adulthood (Fig. 6c and Extended Data Fig. 6c,d). Comparative analyses showed few overlapping rhythmic genes between age groups within each subclass, indicating that the identity of transcripts with rhythmic expression shifts substantially from young adulthood + middle adulthood to late adulthood (Fig. 6c and Supplementary Data 17). Rhythmic genes peaked either shortly after sunrise (at around zeitgeber time 2 (2 ZT); Methods) or in the evening (around 12–14 ZT), a bimodal pattern observed across subclasses and age groups (Fig. 6d). Of the few genes rhythmic in both age groups, 49% (29 out of 59) showed significantly different peak times (FDR < 0.05; Supplementary Data 19), emphasizing age-related reprogramming of rhythmic expression. These findings indicate widespread circadian reprogramming during late adulthood, with both the identity and timing of rhythmic transcripts altered.

In young adulthood + middle adulthood neuronal subclasses, core clock genes displayed synchronized peak expression: forward limb genes peaked at night, while regulatory arms peaked in the morning (Fig. 6e–h, Extended Data Fig. 6b,e, and Supplementary Data 17). In late adulthood, these rhythms were largely lost, and peak times became inconsistent across subclasses (Fig. 6e–h, Extended Data Fig. 6d and Supplementary Data 20). Pathway enrichment analyses mirrored these patterns, with circadian clock pathways highly enriched in young adulthood + middle adulthood but absent in late adulthood (Fig. 6i, Supplementary Fig. 26 and Supplementary Data 21). Notably, microglia and oligodendrocytes in late adulthood showed enrichment for response to unfolded protein, suggesting adaptive transcriptional programs in response to ageing-related stress (Fig. 6i and Supplementary Fig. 26). Further analysis of gene peak timing revealed that sunrise-peaking genes in young adulthood + middle adulthood were enriched for circadian pathways, while sunset-peaking genes in late adulthood were associated with response to unfolded protein in the micro, pericyte (PC), vascular leptomeningeal cell (VLMC) and oligo subclasses (Fig. 6i and Supplementary Data 21). Together, these results demonstrate that, while circadian transcriptional signatures are robust in the young and middle adulthood neuronal subclasses, they are largely lost in late adulthood, with glial subclasses acquiring novel rhythmic programs associated with cellular stress responses during ageing⁴⁴.

Here we present a cell-type-resolved, lifespan-spanning single-nucleus transcriptomic atlas of the human DLPFC, profiling over 1.3 million nuclei from 284 neurotypical donors aged 0 to 97 years. Our atlas reveals that DLPFC nucleus composition follows a non-linear trajectory across the lifespan, with pronounced changes during early development (gliogenesis, synaptic remodelling, apoptosis) followed by stabilization in adulthood. Age 24 emerged as an inflection point, after which most neuronal and glial subclass composition remained stable, except for the IN_SST and OPC subclasses, which continued to show age-associated changes. RNAscope validation confirmed a logarithmic decline in IN_SST and increase in oligodendrocyte abundance across the lifespan.

Our analyses delineate three distinct transcriptomic phases: developmental remodelling, midlife stability and late-life reactivation. During development, 2.9% and 1.7% of aDEGs occurred in neuronal and glial subclasses respectively, stabilizing after age 20 but resurging after age 60 predominantly in glia. While previous studies noted neuronal vulnerability in development⁴⁵ and glial susceptibility with ageing⁴⁶, our findings quantitatively demonstrate that neurons are transcriptionally more resilient to ageing. Ageing modules across glia—microglia, astrocytes, oligodendrocytes—converged on stress-response and immune pathways, suggesting coordinated glial reprogramming in late life.

By quantifying transcriptomic similarity across cell types, we identified developmental divergence reflecting tightly regulated neurodevelopmental processes⁴⁷, contrasting with increasing convergence in late adulthood. Convergence was most pronounced in ENs and enriched for DNA repair and RNA splicing, potentially reflecting shared resilience mechanisms in post-mitotic cells⁴⁸. Thus, we demonstrate subclass-specific transcriptional divergence-to-convergence dynamics across the human DLPFC lifespan.

Modelling continuous gene expression dynamics revealed ten distinct non-linear trajectories spanning all major DLPFC cell types. Trajectory 2, a neuronal resilience program, showed early downregulation of neurodevelopmental disorder risk genes, peaking around ages 12–13 and remaining stable thereafter. This trajectory accounted for 67.6% of developmentally regulated genes and aligns with normative brain growth patterns, suggesting that perturbations before this critical window may have long-term functional impacts⁴⁹. By contrast, trajectory 10 captured a glial ageing program defined by late-life upregulation of immune and unfolded protein response genes.

Pseudotime analyses identified dynamic transcriptional programs and disease associations. Psychiatric-trait-linked genes showed persistent neuronal expression across the lifespan, indicating roles in both development and adult brain maintenance⁴⁵, with high expression during astrocyte and oligodendrocyte development. Convergence in ENs during late adulthood associates with DNA repair and RNA splicing to preserve post-mitotic cell integrity. Conversely, neurodegenerative disease-associated genes were highly expressed in ageing microglia and oligodendrocytes, suggesting a shift in glial functions from developmental support to late-life homeostasis and immunity^10,34. These late-life glial programs implicate shared ageing mechanisms underlying diverse neurodegenerative conditions.

Spatial transcriptomics validation showed EN modules localized to expected cortical layers, while glial modules exhibited grey-to-white matter enrichment during ageing, supporting predicted compartmental transitions and highlighting region-specific vulnerability. These integrated analyses provide an anatomically resolved, lifespan-wide map of human DLPFC gene regulatory programs.

Circadian rhythm analysis identified marked reprogramming in late adulthood. Although core clock gene rhythmicity was robust in young and middle adulthood, it diminished in neurons during late adulthood. Conversely, microglia and oligodendrocytes gained rhythmicity in unfolded protein response pathways, suggesting adaptive stress responses^42,43. We present a single-cell characterization of circadian gene expression across adult stages of the human DLPFC.

Our study is limited by reliance on transcriptomic data alone, which does not capture post-transcriptional, proteomic or epigenomic regulation. Future integration with proteomics, epigenomics and multiregion spatial transcriptomics will elucidate how transcriptional ageing trajectories translate to functional changes and clarify how resilience mechanisms in development relate to vulnerability in ageing.

This atlas provides a foundation for understanding cellular program transitions from resilience to vulnerability. Our findings highlight neuronal transcriptional stability across adulthood, late-life glial reprogramming, and the interplay between ageing, circadian disruption and disease risk. Distinct developmental and ageing inflection points represent windows for therapeutic intervention, offering a critical resource for targeted strategies to preserve brain health and mitigate age-related cognitive decline.

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 <cut-off and ≥cut-off groups. This revealed a plateau in age-related changes after the age of 23 years, suggesting stabilization in adulthood. We then selected the age of 24 years as an optimal breakpoint and visualized subclass-specific trends in development versus adulthood.

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.

All the source codes used in this study are available at GitHub (https://github.com/DiseaseNeuroGenomics/DLPFC_Lifespan_Atlas).

1. Girdhar, K. et al. Chromatin domain alterations linked to 3D genome organization in a large cohort of schizophrenia and bipolar disorder brains. Nat. Neurosci. 25, 474–483 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

2. Cain, A. et al. Multicellular communities are perturbed in the aging human brain and Alzheimer’s disease. Nat. Neurosci. 26, 1267–1280 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

3. Nejati, V., Majdi, R., Salehinejad, M. A. & Nitsche, M. A. The role of dorsolateral and ventromedial prefrontal cortex in the processing of emotional dimensions. Sci. Rep. 11, 1971 (2021).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

4. Velmeshev, D. et al. Single-cell analysis of prenatal and postnatal human cortical development. Science 382, eadf0834 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

5. Zeng, B. et al. The single-cell and spatial transcriptional landscape of human gastrulation and early brain development. Cell Stem Cell 30, 851–866 (2023).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

6. Herring, C. A. et al. Human prefrontal cortex gene regulatory dynamics from gestation to adulthood at single-cell resolution. Cell 185, 4428–4447 (2022).

   Article  CAS  PubMed  Google Scholar 

7. Emani, P. S. et al. Single-cell genomics and regulatory networks for 388 human brains. Science 384, eadi5199 (2024).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

8. Fullard, J. F. et al. Population-scale cross-disorder atlas of the human prefrontal cortex at single-cell resolution. Sci. Data https://doi.org/10.1038/s41597-025-04687-5 (2025).

9. Lee, D. et al. Single-cell atlas of transcriptomic vulnerability across brain disorders. Nature https://doi.org/10.1038/s41586-025-09573-z (2026).

10. Mathys, H. et al. Single-cell atlas reveals correlates of high cognitive function, dementia, and resilience to Alzheimer’s disease pathology. Cell 186, 4365–4385 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

11. Hoffman, G. E. & Schadt, E. E. variancePartition: interpreting drivers of variation in complex gene expression studies. BMC Bioinform. 17, 483 (2016).

    Article  Google Scholar 

12. Chaudhari, P. R., Singla, A. & Vaidya, V. A. Early adversity and accelerated brain aging: a mini-review. Front. Mol. Neurosci. 15, 822917 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

13. Felsky, D. et al. Polygenic analysis of inflammatory disease variants and effects on microglia in the aging brain. Mol. Neurodegener. 13, 38 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

14. Leutner, M. et al. Obesity as pleiotropic risk state for metabolic and mental health throughout life. Transl. Psychiatry 13, 175 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

15. Hoffman, G. E. & Roussos, P. Fast, flexible analysis of differences in cellular composition with crumblr. Nat. Commun. https://doi.org/10.1038/s41467-026-75681-7 (2026).

16. Hoffman, G. E. et al. Efficient differential expression analysis of large-scale single-cell transcriptomics data using Dreamlet. Nat. Commun. https://doi.org/10.1038/s41467-026-75680-8 (2026).

17. Ling, E. et al. A concerted neuron-astrocyte program declines in ageing and schizophrenia. Nature 627, 604–611 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

18. Ruzicka, W. B. et al. Single-cell multi-cohort dissection of the schizophrenia transcriptome. Science 384, eadg5136 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

19. Karczewski, K. J. et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature 581, 434–443 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

20. Lynch, G., Rex, C. S. & Gall, C. M. Synaptic plasticity in early aging. Ageing Res. Rev. 5, 255–280 (2006).

    Article  CAS  PubMed  Google Scholar 

21. Arcos-Burgos, M., Lopera, F., Sepulveda-Falla, D. & Mastronardi, C. Neural plasticity during aging. Neural Plast. 2019, 6042132 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

22. Izgi, H. et al. Inter-tissue convergence of gene expression during ageing suggests age-related loss of tissue and cellular identity. eLife 11, e68048 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Jin, K. et al. Brain-wide cell-type-specific transcriptomic signatures of healthy ageing in mice. Nature 638, 182–196 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

24. Urbut, S. M., Wang, G., Carbonetto, P. & Stephens, M. Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nat. Genet. 51, 187–195 (2019).

    Article  CAS  PubMed  Google Scholar 

25. Lüleci, H. B. & Yılmaz, A. Robust and rigorous identification of tissue-specific genes by statistically extending tau score. BioData Min. 15, 31 (2022).

    Article  PubMed  PubMed Central  Google Scholar 

26. Köhler, S., Winkler, U. & Hirrlinger, J. Heterogeneity of astrocytes in grey and white matter. Neurochem. Res. 46, 3–14 (2021).

    Article  PubMed  Google Scholar 

27. Kruyer, A., Kalivas, P. W. & Scofield, M. D. Astrocyte regulation of synaptic signaling in psychiatric disorders. Neuropsychopharmacology 48, 21–36 (2023).

    Article  PubMed  Google Scholar 

28. Linnerbauer, M., Wheeler, M. A. & Quintana, F. J. Astrocyte crosstalk in CNS inflammation. Neuron 108, 608–622 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Romanos, J., Benke, D., Pietrobon, D., Zeilhofer, H. U. & Santello, M. Astrocyte dysfunction increases cortical dendritic excitability and promotes cranial pain in familial migraine. Sci. Adv. 6, eaaz1584 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

30. Sofroniew, M. V. & Vinters, H. V. Astrocytes: biology and pathology. Acta Neuropathol. 119, 7–35 (2010).

    Article  PubMed  Google Scholar 

31. Badia-I-Mompel, P. et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform. Adv. 2, vbac016 (2022).

    Article  PubMed  PubMed Central  Google Scholar 

32. Wang, L. et al. Molecular and cellular dynamics of the developing human neocortex. Nature 1, 10 (2025).

    Google Scholar 

33. Huuki-Myers, L. A. et al. A data-driven single-cell and spatial transcriptomic map of the human prefrontal cortex. Science 384, eadh1938 (2024).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

34. Bartels, T., De Schepper, S. & Hong, S. Microglia modulate neurodegeneration in Alzheimer’s and Parkinson’s diseases. Science 370, 66–69 (2020).

    Article  ADS  CAS  PubMed  Google Scholar 

35. Michell-Robinson, M. A. et al. Roles of microglia in brain development, tissue maintenance and repair. Brain 138, 1138–1159 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

36. Bellenguez, C. et al. New insights into the genetic etiology of Alzheimer’s disease and related dementias. Nat. Genet. 54, 412–436 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

37. Greig, L. C., Woodworth, M. B., Galazo, M. J., Padmanabhan, H. & Macklis, J. D. Molecular logic of neocortical projection neuron specification, development and diversity. Nat. Rev. Neurosci. 14, 755–769 (2013).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

38. Südhof, T. C. The cell biology of synapse formation. J. Cell Biol. 220, e202103052 (2021).

39. Forsyth, J. K. & Lewis, D. A. Mapping the consequences of impaired synaptic plasticity in schizophrenia through development: an integrative model for diverse clinical features. Trends Cogn. Sci. 21, 760–778 (2017).

    Article  PubMed  PubMed Central  Google Scholar 

40. Johnson, K. A. et al. Tourette syndrome: clinical features, pathophysiology, and treatment. Lancet Neurol. 22, 147–158 (2023).

    Article  PubMed  Google Scholar 

41. Lowe, C. J., Reichelt, A. C. & Hall, P. A. The prefrontal cortex and obesity: a health neuroscience perspective. Trends Cogn. Sci. 23, 349–361 (2019).

    Article  PubMed  Google Scholar 

42. Stowe, T. A. & McClung, C. A. How does chronobiology contribute to the development of diseases in later life. Clin. Interv. Aging 18, 655–666 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

43. Wolff, C. A. et al. Defining the age-dependent and tissue-specific circadian transcriptome in male mice. Cell Rep. 42, 111982 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

44. Uddin, M. S., Yu, W. S. & Lim, L. W. Exploring ER stress response in cellular aging and neuroinflammation in Alzheimer’s disease. Ageing Res. Rev. 70, 101417 (2021).

    Article  CAS  PubMed  Google Scholar 

45. Li, M. et al. Integrative functional genomic analysis of human brain development and neuropsychiatric risks. Science 362, eaat7615 (2018).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

46. Green, G. S. et al. Cellular dynamics across aged human brains uncover a multicellular cascade leading to Alzheimer’s disease. Preprint at bioRxiv https://doi.org/10.1101/2023.03.07.531493 (2023).

47. Nowakowski, T. J. et al. Spatiotemporal gene expression trajectories reveal developmental hierarchies of the human cortex. Science 358, 1318–1323 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

48. Reid, D. A. et al. Incorporation of a nucleoside analog maps genome repair sites in postmitotic human neurons. Science 372, 91–94 (2021).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

49. Khodosevich, K. & Sellgren, C. M. Neurodevelopmental disorders-high-resolution rethinking of disease modeling. Mol. Psychiatry 28, 34–43 (2023).

    Article  PubMed  Google Scholar 

50. Stoeckius, M. et al. Cell Hashing with barcoded antibodies enables multiplexing and doublet detection for single cell genomics. Genome Biol. 19, 224 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

51. Kaminow, B., Yunusov, D. & Dobin, A. STARsolo: accurate, fast and versatile mapping/quantification of single-cell and single-nucleus RNA-seq data. Preprint at bioRxiv https://doi.org/10.1101/2021.05.05.442755 (2021).

52. Huang, X. & Huang, Y. Cellsnp-lite: an efficient tool for genotyping single cells. Bioinformatics https://doi.org/10.1093/bioinformatics/btab358 (2021).

53. Huang, Y., McCarthy, D. J. & Stegle, O. Vireo: Bayesian demultiplexing of pooled single-cell RNA-seq data without genotype reference. Genome Biol. 20, 273 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

54. Fort, A. et al. MBV: a method to solve sample mislabeling and detect technical bias in large combined genotype and sequencing assay datasets. Bioinformatics 33, 1895–1897 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

55. Li, B. et al. Cumulus provides cloud-based data analysis for large-scale single-cell and single-nucleus RNA-seq. Nat. Methods 17, 793–798 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

56. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

57. Fleming, S. J. et al. Unsupervised removal of systematic background noise from droplet-based single-cell experiments using CellBender. Nat. Methods 20, 1323–1335 (2023).

    Article  CAS  PubMed  Google Scholar 

58. Wolock, S. L., Lopez, R. & Klein, A. M. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 8, 281–291 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

59. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 16, 1289–1296 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Xu, C. et al. Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models. Mol. Syst. Biol. 17, e9620 (2021).

    Article  PubMed  PubMed Central  Google Scholar 

61. de Leeuw, C. A., Mooij, J. M., Heskes, T. & Posthuma, D. MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput. Biol. 11, e1004219 (2015).

    Article  PubMed  PubMed Central  Google Scholar 

62. Gayoso, A. et al. A Python library for probabilistic analysis of single-cell omics data. Nat. Biotechnol. 40, 163–166 (2022).

    Article  CAS  PubMed  Google Scholar 

63. Szklarczyk, D. et al. The STRING database in 2023: protein–protein association networks and functional enrichment analyses for any sequenced genome of interest. Nucleic Acids Res. 51, D638–D646 (2022).

    Article  Google Scholar 

64. Pardo, B. et al. spatialLIBD: an R/Bioconductor package to visualize spatially-resolved transcriptomics data. BMC Genom. 23, 434 (2022).

    Article  CAS  Google Scholar 

65. 1000 Genomes Project Consortium. A global reference for human genetic variation. Nature 526, 68–74 (2015).

    Article  Google Scholar 

66. Müller-Dott, S. et al. Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Res. 51, 10934–10949 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

67. Seney, M. L. et al. Diurnal rhythms in gene expression in the prefrontal cortex in schizophrenia. Nat. Commun. 10, 3355 (2019).

    Article  ADS  PubMed  PubMed Central  Google Scholar 

68. Xue, X. et al. DiffCircaPipeline: a framework for multifaceted characterization of differential rhythmicity. Bioinformatics 39, btad039 (2023).

69. Kim, Y. H. & Lazar, M. A. Transcriptional control of circadian rhythms and metabolism: a matter of time and space. Endocr. Rev. 41, 707–732 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Author notes

1. These authors contributed equally: Hui Yang, Tereza Clarence, Madeline R. Scott, Kiran Girdhar, Panos Roussos

Authors and Affiliations

1. Center for Disease Neurogenomics, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Hui Yang, Tereza Clarence, Xinyi Wang, Prashant N. M., Milos Pjanic, Sanan Venkatesh, Aram Hong, Clara Casey, Sarah R. Murphy, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Alexander Kawah Yu, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, David Burstein, Deepika Mathur, Fotios Tsetsos, Jennifer Monteiro Fortes, Karen Therrien, Lyra Sheu, Marios Anyfantakis, Mikaela Koutrouli, Nicolas Y. Masse, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Seon Kinrot, Steven P. Kleopoulos, Zhenyi Wu, Georgios Voloudakis, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Kiran Girdhar & Panos Roussos

2. Friedman Brain Institute, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Hui Yang, Tereza Clarence, Xinyi Wang, Prashant N. M., Milos Pjanic, Sanan Venkatesh, Aram Hong, Clara Casey, Sarah R. Murphy, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Alexander Kawah Yu, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, David Burstein, Deepika Mathur, Fotios Tsetsos, Jennifer Monteiro Fortes, Karen Therrien, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Mikaela Koutrouli, Nicolas Y. Masse, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Seon Kinrot, Steven P. Kleopoulos, Zhenyi Wu, Nadejda M. Tsankova, Vahram Haroutunian, Georgios Voloudakis, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Kiran Girdhar & Panos Roussos

3. Department of Psychiatry, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Hui Yang, Tereza Clarence, Xinyi Wang, Prashant N. M., Milos Pjanic, Sanan Venkatesh, Aram Hong, Clara Casey, Sarah R. Murphy, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Alexander Kawah Yu, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, David Burstein, Deepika Mathur, Fotios Tsetsos, Jennifer Monteiro Fortes, Karen Therrien, Lyra Sheu, Marios Anyfantakis, Maxim Signaevsky, Mikaela Koutrouli, Nicolas Y. Masse, Pavel Katsel, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Seon Kinrot, Steven P. Kleopoulos, Zhenyi Wu, Vahram Haroutunian, Georgios Voloudakis, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Kiran Girdhar & Panos Roussos

4. Department of Genetics and Genomic Science, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Hui Yang, Tereza Clarence, Xinyi Wang, Prashant N. M., Milos Pjanic, Sanan Venkatesh, Aram Hong, Clara Casey, Sarah R. Murphy, Zhiping Shao, Marcela Alvia, Stathis Argyriou, Alexander Kawah Yu, Biao Zeng, Christian Dillard, Christian Porras, Collin Spencer, David Burstein, Deepika Mathur, Fotios Tsetsos, Jennifer Monteiro Fortes, Karen Therrien, Lyra Sheu, Marios Anyfantakis, Mikaela Koutrouli, Nicolas Y. Masse, Pengfei Dong, Rachel Bercovitch, Roman Kosoy, Seon Kinrot, Steven P. Kleopoulos, Zhenyi Wu, Georgios Voloudakis, Jaroslav Bendl, Donghoon Lee, John F. Fullard, Gabriel E. Hoffman, Kiran Girdhar & Panos Roussos

5. Department of Psychiatry, University of Pittsburgh School of Medicine, Pittsburgh, PA, USA

   Madeline R. Scott & Colleen A. McClung

6. Department of Neuroscience, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Nadejda M. Tsankova & Vahram Haroutunian

7. Department of Pathology, Molecular, and Cell-Based Medicine, Icahn School of Medicine at Mount Sinai, New York, NY, USA

   Nadejda M. Tsankova

8. Human Brain Collection Core, National Institute of Mental Health-Intramural Research Program, Bethesda, MD, USA

   Pavan K. Auluck & Stefano Marenco

9. Mental Illness Research, Education and Clinical Center VISN2, James J. Peters VA Medical Center, New York, NY, USA

   David Burstein, Deepika Mathur, Vahram Haroutunian, Georgios Voloudakis, Gabriel E. Hoffman & Panos Roussos

10. Department of Artificial Intelligence and Human Health, Icahn School of Medicine at Mount Sinai, New York, NY, USA

    Georgios Voloudakis

11. Center for Precision Medicine and Translational Therapeutics, James J. Peters VA Medical Center, New York, NY, USA

    Georgios Voloudakis, Gabriel E. Hoffman & Panos Roussos

12. Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA

    Athan Z. Li, Daifeng Wang, Noah Cohen Kalafut & Sayali A. Alatkar

13. Waisman Center, University of Wisconsin-Madison, Madison, WI, USA

    Athan Z. Li, Chenfeng He, Chirag Gupta, Daifeng Wang, Jerome J. Choi, Kalpana H. Arachchilage, Noah Cohen Kalafut, Pramod B. Chandrashekar, Saniya Khullar, Sayali A. Alatkar, Ting Jin & Xiang Huang

14. Department of Biostatistics and Medical Informatics, University of Wisconsin-Madison, Madison, WI, USA

    Chenfeng He, Chirag Gupta, Daifeng Wang, Kalpana H. Arachchilage, Pramod B. Chandrashekar, Saniya Khullar & Ting Jin

15. Rush Alzheimer’s Disease Center, Rush University Medical Center, Chicago, IL, USA

    David A. Bennett, David Burstein, Deepika Mathur, Fotios Tsetsos & Lisa L. Barnes

16. Department of Neurological Sciences, Rush University Medical Center, Chicago, IL, USA

    David A. Bennett & Lisa L. Barnes

17. Center for Systems and Therapeutics, Gladstone Institutes, San Francisco, CA, USA

    Gennadi Ryan, Monika Ahirwar & Vivek G. Ramaswamy

18. Department of Physiology, University of California San Francisco, San Francisco, CA, USA

    Gennadi Ryan & Steven Finkbeiner

19. Department of Population Health Sciences, University of Wisconsin-Madison, Madison, WI, USA

    Jerome J. Choi & Steven Finkbeiner

20. Novo Nordisk Foundation Center for Protein Research, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark

    Lars J. Jensen & Mikaela Koutrouli

21. Vanderbilt Genetics Institute, Vanderbilt University Medical Center, Nashville, TN, USA

    Logan C. Dumitrescu & Timothy J. Hohman

22. Vanderbilt Memory & Alzheimer’s Center, Vanderbilt University Medical Center, Nashville, TN, USA

    Logan C. Dumitrescu & Timothy J. Hohman

23. Taube/Koret Center for Neurodegenerative Disease Research, Gladstone Institutes, San Francisco, CA, USA

    Monika Ahirwar, Steven Finkbeiner & Vivek G. Ramaswamy

24. Neuroscience and Biomedical Sciences Graduate Programs, University of California San Francisco, San Francisco, CA, USA

    Steven Finkbeiner

25. Department of Neurology, University of California San Francisco, San Francisco, CA, USA

    Steven Finkbeiner

Consortia

Contributions

Conceptualization: K.G. and P.R. Methodology and software: G.E.H., K.G., H.Y. and M.R.S. Validation: H.Y., K.G. and X.W. Formal analysis: K.G., H.Y., M.R.S. and T.C. Investigation: J.F.F., X.W., M.P., S.V., A.H., C.C., S.R.M., Z.S., M.A., S.A. and A.K.Y. Resources: V.H., P.K.A., S.M. and N.M.T. Data collection: P.R., V.H., S.F. and D.W. Data curation: J.B., P.N.M., S.V., M.P., D.L. and H.Y. Writing: H.Y., K.G., M.R.S., T.C., P.R., J.F.F., X.W. and D.L. with support from all of the other authors. Visualization: K.G., H.Y., M.R.S. and T.C. Supervision: P.R., K.G., J.F.F., C.A.M., G.V., J.B., V.H., D.L. and G.E.H. Project administration: K.G. and P.R. Funding acquisition: P.R. and V.H. All of the authors read and approved the final draft of the paper.

Corresponding authors

Correspondence to Kiran Girdhar or Panos Roussos.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Qin Ma and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

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