Human brain organoids record the passage of time over multiple years

Human brain development is orchestrated by transcriptional programs defined by spatially and temporally regulated waves of gene expression. These gene expression patterns are tightly coordinated by dynamic changes in the activity of regulatory elements and epigenomic remodelling, ensuring the timely emergence of distinct neural cell types and the progressive maturation of the nervous system. In humans, brain development and maturation proceed at a much slower pace than in most other species, with cortical neurons requiring years to reach full maturity. This protracted timeline is preserved in all cell types generated in vitro from human pluripotent stem (hPS) cells, including cortical neurons, pointing at a cell-intrinsic clock that sets the pace of brain development, although its molecular mechanisms and functional significance remain to be elucidated^1,2. Recent work has shown the importance of epigenetic barriers that enforce the slow timing of human neuronal maturation¹, and the role of species-specific rates of mitochondrial metabolism in influencing developmental tempo³. Human organoids could, in principle, be used to study these processes of human brain maturation, but existing models largely recapitulate earlier developmental processes, and most studies do not follow organoids over a sufficiently long time span. One report cultured organoids for up to 694 days, and used bulk RNA-sequencing (RNA-seq) and methylation arrays to profile whole organoids⁴. However, there is known variation in the ability of distinct cell types to survive over long periods in culture, especially neuronal populations. Moreover, beyond the continued presence of a cell type, there is a need to understand whether structural and functional properties such as neuronal architecture and coordinated circuit activity are present over years in culture. It is therefore important to resolve and analyse individual cell populations and their functional features across different modalities along these extended timelines. Without cell-type-specific knowledge, and in the absence of longer timelines of maturation, we currently lack understanding of the mechanisms by which the many individual cell types of the human brain mature in vitro, and whether they can measure and record time as they do in vivo.

Here we developed human cortical organoids for over 5 years in culture and integrated single-cell transcriptional information with epigenetic, structural and functional data to build a comprehensive map of development and maturation across an unprecedented time span. The data indicate that human brain organoids record the passage of time following endogenous milestones and using similar epigenetic mechanisms. Importantly, progenitors in old organoids can recall the passage of time, like a memory of development that has already happened, to generate late progeny without repeating early developmental steps. The work demonstrates that a diversity of cells in organoids can develop over considerable time periods and are capable of both recording and recalling their developmental age. These systems provide a wealth of data on largely inscrutable periods of postnatal development of the human brain.

The development and maturation of human brain organoids follows the notoriously slow, neotenous pace of the endogenous human brain, which requires many years to reach adulthood^5,6,7,8. To date, organoids have mostly been studied over timespans of a few months, and it remains unclear whether they can continue to mature and age in culture over many years to mimic aspects of postnatal maturation of the human brain that otherwise remain largely experimentally inaccessible. To investigate this possibility, we cultured human organoids of the cerebral cortex for over 5 years using a protocol that we published previously⁷. To characterize the cellular composition and molecular features of individual cell types in these long-term cultures, we used single-cell RNA-seq (scRNA-seq) analysis to profile 34 single organoids collected at 6 months, 9 months, and 1, 1.5, 2, 3, 4 and 5 years, and integrated these data with our previously published dataset of 76 single organoids profiled from 15 days to 6 months^6,7,9 (Fig. 1a,b, Methods and Supplementary Table 1; n = 110 organoids in total, from 34 new datasets and 76 previously generated datasets; n = 424,720 cells).

To correlate transcriptional changes in organoids over years in vitro with those occurring in the human brain over years in vivo, we used label transfer from reference datasets of endogenous human cortex^10,11 to estimate the transcriptional ‘age’ of organoid cells (Methods): early organoids (15 days to 2 months) mapped to first-trimester human fetal brains; 3 to 6 month organoids predominantly mapped to the second trimester; and the later stages (9 months to 5 years) progressively shifted towards mapping to late prenatal and postnatal ages (Fig. 1c). This suggested a temporal alignment with in vivo developmental trajectories, with postnatal-like signatures emerging in organoids after 12 months in culture.

To precisely assess transcriptional maturation, we used DIALOGUE¹² to identify gene modules correlated with age in endogenous human fetal and postnatal cortex datasets^10,11 (Fig. 1d–f, Extended Data Fig. 1a,b, Supplementary Information and Supplementary Table 2). Projecting these maturation modules onto the organoid time-series dataset revealed a strong correlation between in vivo age-related module scores and duration of in vitro development (Pearson’s r = −0.08, 0.78, 0.35, 0.78 and −0.72 for multicellular programs (MCPs) 1–5, respectively), showing that cells in organoids preserve a sequential progression of gene expression programs associated with maturation in vivo.

DNA methylation, particularly at age-associated CpG sites, serves as a powerful molecular readout of biological age¹³. In vivo, in cortex and other tissues, the methylome is dynamic during early development, but eventually reaches high levels of DNA methylation that are stably maintained postnatally. Whole-genome bisulfite sequencing (WGBS) of single organoids at nine timepoints from 3 months to 5 years in culture showed that organoids follow a similar global pattern of progressive methylation and subsequent stability, with dynamics and features that resemble epigenetic features of endogenous tissues. Key cortical genes involved in fate acquisition and layer specification showed progressive, time-dependent gains in CpG methylation, particularly within CpG islands and DNA methylation valleys (DMVs; Fig. 1g,h, Extended Data Fig. 1c–g and Supplementary Tables 3 and 4; a detailed discussion is provided in the Supplementary Information).

When comparing with endogenous postnatal maturation programs, we found that cell-type-specific differentially methylated regions (cdDMRs) that are variable in neuronal and non-neuronal cells across human cortical development from 0 to 23 years¹⁴ likewise vary in organoids, in a manner that is largely consistent with their in vivo dynamics (Extended Data Fig. 1h–k and Supplementary Information). Differential methylation analysis between the 3-month and 5-years timepoints identified 213 age-correlated DMRs, which were enriched in epigenetic features of known developmental importance, including CpG islands, DMVs and gene regulatory elements, and included many genes known to have key roles in the developing neocortex (Fig. 1i, Extended Data Fig. 2a–d, Supplementary Table 3 and Supplementary Information). Finally, organoids showed progressive CpA methylation (mCA; a non-CpG modification that is a hallmark of neuronal maturity¹⁵) (Extended Data Fig. 2e) at a set of 1,383 human fetal forebrain superenhancers that gain mCA during postnatal neuron development in vivo¹⁵, indicating appropriate temporal regulation of mCA programs (Extended Data Fig. 2f,g).

Together, these results imply that the ongoing methylation dynamics observed in organoids over extended culture may reflect in vitro recapitulation of endogenous development and maturation epigenetic programs.

CpG sites that show reproducible gains or losses in methylation over time¹⁶ have been used to create age clocks, correlating methylation with the temporal age of the tissue, in vivo¹⁷. To evaluate whether organoids capture features of biological ageing, we applied three epigenetic clocks: the Horvath pan-tissue clock¹⁷, a human fetal brain-specific epigenetic clock¹⁸ and a human cortex-specific DNA methylation clock (Fig. 1j–k). In both the Horvath and cortex-specific models, predicted DNA methylation age (DNAm age) tracked closely with culture time (r = 0.88–0.90; median absolute error: 7.25 and 20.04 months, respectively) (Fig. 1j,k). Similarly, estimates from the fetal brain clock showed positive relationship with culture duration (r = 0.54, P = 0.056; Extended Data Fig. 2h). To verify that the observed methylation changes were not primarily driven by proliferation, we confirmed that methylation at solo-WCGW tetranucleotide sites, which correlates with mitotic division history, remained stable (Extended Data Fig. 2i).

Collectively, these findings demonstrate that organoids are capable of preserving robust, time-resolved transcriptional and DNA methylation programs that are associated with tissue ageing in vivo. The ability to recapitulate transcriptional and epigenetic signatures of in vivo maturation across 5 years in vitro establishes organoids as a powerful model system to study the biology of both the prenatal and postnatal human brain.

It is well known that not all cell types survive equally well in organoids, underscoring the importance of single-cell resolution when assaying long-term organoid cultures. Our scRNA-seq time course showed that all expected cell types were generated in cortical organoids (Fig. 2a,b, Extended Data Fig. 3a,b and Supplementary Table 5). Most broad cell types exhibited a significant correlation between age and MCP4 values, indicating that each population continued to mature over time in culture (Extended Data Fig. 3d,e). Reproducibility metrics, assessed using Aitchison distances, remained within the range of human tissue variability throughout the entire culture period (Extended Data Fig. 3c).

We next investigated whether individual cell types in organoids maintained molecular fidelity to the human reference over time. Across astrocytes and glial progenitors, which remained well-represented over the 5 years of culture, we found consistent mapping of time in vitro to progressively later human developmental ages (Extended Data Fig. 3f). Of the genes that showed statistically significant change over time in both endogenous and organoid glial populations (1,907 genes), the majority (80.2%) displayed concordant temporal trends between organoids and endogenous tissue, and no significantly enriched biological processes were detected in the minority (19.8%) that showed discordant correlation (Extended Data Fig. 3g). Importantly, gene set enrichment analysis (GSEA) revealed that the glial transcriptional changes in vivo and in vitro shared underlying molecular signatures (Extended Data Fig. 3h). Finally, DIALOGUE revealed the involvement of cell-type-specific genes with known developmental roles in age-correlated multicellular processes (Extended Data Fig. 4a). Consistent with this molecular evidence of glial maturation, organoids gained progressive evidence of myelination at both the ultrastructural and protein levels (Extended Data Fig. 4b,c), indicating the presence of myelinating oligodendrocytes by 6 months in culture.

We next considered neurons. In this culture condition (CDM4; Methods), we find that neuronal populations progressively declined in proportion over time in the scRNA-seq data, although we still detected a small cluster containing inhibitory and excitatory neuronal markers in the 5-year-old samples (Fig. 2a,b and Extended Data Fig. 3a,b). As the dissociation required for scRNA-seq may lead to differential death of fragile cells such as neurons, we assessed the presence of neurons using methods applicable to intact organoids. Cells expressing the neuronal marker NeuN were detected by immunohistochemistry at 2, 3, 4 and 5.8 years in culture, and the excitatory neuronal marker SATB2 at 5.8 years. In agreement, bulk RNA-seq analysis of whole organoids demonstrated expression of neuronal genes at 2 and 3 years (Extended Data Fig. 2j,k). Together, these results demonstrate that despite a progressive decline in neuronal representation over time, long-term organoids retain surviving neuronal populations even after nearly 6 years in culture. We conclude that while these culture conditions can maintain neurons over a span of years, they would benefit from further optimization to ensure consistent representation of neurons across extended timelines.

Spontaneous firing activity is known to have a role in neuronal survival, but conventional medium does not support robust activity. We therefore investigated a previously established medium designed to preserve spontaneous activity (BrainPhys)¹⁹, to test whether it could improve neuronal survival in organoids. We used BrainPhys medium supplemented with GlutaMax to improve stability of glutamine and reduce ammonia accumulation²⁰ (hereafter, activity permissive medium (APM); Methods), and cultured organoids in APM starting from day in vitro 70 (DIV70; Fig. 2c). At 9 months, organoids cultured in APM showed a significant increase in the number of both callosal projection neurons (CPNs) (SATB2⁺ cells) (n = 9 APM and 9 CDM4 organoids; linear mixed model (LMM), P = 2.4 × 10⁻⁵, 95% confidence interval (CI) = 68.1–144.4) and active CPNs (SATB2⁺FOS⁺ double-positive cells) (n = 9 APM and 9 CDM4 organoids; LMM, P = 9.0 × 10⁻⁴, 95% CI = 16.6–52.3) (Fig. 2d), consistent with an increased number of active neurons in the APM organoids.

To investigate whether APM promoted ultrastructural and morphological features associated with neuronal maturity and activity, we examined synapse density and synaptic spine formation using quantitative electron microscopy (EM) and neuronal morphology by expansion microscopy. EM with automated synapse detection²¹ showed that, while no statistically significant differences in synaptic density were observed at 6 months (Wilcoxon rank-sum test, P = 0.7), at 1 year, APM organoids exhibited higher synaptic densities (P = 0.002) than CDM4 organoids (6 months: n = 3 organoids from each of n = 2 genetic backgrounds (H1 and 11a); 1 year: n = 6 organoids from each of n = 3 genetic backgrounds (H1, 11a and PGP1), for both CDM4 and APM conditions; Fig. 2e–g, Supplementary Table 6, Supplementary Videos 1 and 2 and Supplementary Information). Furthermore, at 1 year, a significantly higher proportion of synapses in the APM-treated organoids was located on spines (25% versus 52% in CDM4 and APM, respectively; two-sided Fisher’s exact tests: 6 months: P = 0.51, odds ratio (OR) = 0.88, 95% CI = 0.60–1.30; 1 year: P = 9.5 × 10⁻⁸, OR = 4.67, 95% CI = 2.49–9.23; Fig. 2f and Supplementary Information). This is consistent with enhanced synaptic maturation and stabilization in APM organoids, features that are critical for the long-term integrity of synaptic circuits.

To assess neuronal complexity (that is, arborization and process outgrowth)—a clear measure of cortical neuron maturation in vivo—we used epitope tag barcoding^22,23,24,25 as a tracing strategy²⁶ in combination with antigenically orthogonal fluorescent proteins²⁷ and magnify expansion microscopy²⁸ (Fig. 2h,i, Extended Data Fig. 5a,b and Methods). At 7 months, APM-treated neurons exhibited significantly greater morphological complexity (total neurite length, maximal process extension and number of terminal branches), along with significantly higher maximum Strahler numbers—a measure of arborization complexity (Fig. 2j–n; a detailed discussion is provided in the Supplementary Information). Together, these results demonstrate that APM culture conditions promote both synapse development and increased neuronal morphological complexity.

To examine the effects of APM on cell proportions and transcriptional states over time, we performed scRNA-seq analysis of individual organoids from four different donors over a time course from 4 months to 1.5 years in culture (4 months: CDM4: n = 6, 15,972 cells; APM: n = 2, 4,579 cells; 6 months: CDM4: n = 10, 64,817 cells; APM: n = 9, 41,231 cells; 9 months: CDM4: n = 6, 18,104 cells; APM: n = 3, 8,533 cells; 1 year: CDM4: n = 4, 7,802 cells; APM: n = 3, 6,084 cells; 1.5 years: CDM4: n = 2, 3,952 cells; APM: n = 2, 8,670 cells; Fig. 3a–g and Extended Data Figs. 5c–e and 6 and Supplementary Table 7). Reproducibility, assessed using Atchison distances⁹, demonstrated that, at all ages in both conditions, intersample variability fell within the limits of variation found in endogenous tissue (Fig. 3f and Extended Data Fig. 6j).

We examined the maintenance of excitatory neuronal populations, with particular focus on CPNs, a late-born population of excitatory neurons that are expanded in humans compared with other mammals and are associated with the evolution of higher cognitive functions²⁹. Marked differences in cellular composition were apparent as early as 6 months (CDM4: n = 10 organoids, 64,817 cells; APM: n = 9 organoids, 41,231 cells): APM organoids showed increased proportions of excitatory neuron subtypes, including CPNs and corticofugal projection neurons (CFuPNs), while glial populations and inhibitory neurons were proportionately reduced (false-discovery rate (FDR)-adjusted P = 1.15 × 10⁻² (apical radial glia, aRGs), 3.38 × 10⁻⁷ (outer radial glia, oRGs), 1.95 × 10⁻² (intermediate progenitors, IPs), 9.23 × 10⁻⁸ (CFuPNs), 2.69 × 10⁻⁸ (CPNs), 1.24 × 10⁻⁵ (interneuron (IN) progenitors), 4.19 × 10⁻⁶ (immature INs), 4.81 × 10⁻² (truncated radial glia, tRGs), 1.52 × 10⁻¹ (glial precursors), 4.33 × 10⁻¹³ (astrocytes), 2.46 × 10⁻² (oligodendrocyte progenitor cells, OPCs), negative binomial mixed-effect (NBME) model; Fig. 3b–g). Notably, the enrichment of CPNs was maintained at 9 months and 1 year (FDR-adjusted P = 7.05 × 10⁻³ (9 months) and 1.96 × 10⁻² (1 year)) and, at 1.5 years, we observed an absolute increase in CPNs (3.04%) in APM compared with CDM4 (Fig. 3b and Extended Data Fig. 6). Together, this indicates that APM culture conditions selectively favour the long-term maintenance and stability of excitatory cortical neurons, including CPNs.

Consistent with our findings in EM and morphological analysis, neurons from APM organoids showed increased gene module scores for synapse-related genes (6 months: n = 9 APM and 10 CDM4 organoids; LMM, P = 0.022, expected increase = 1.03, 95% CI = 1.01–1.05; Methods, Extended Data Fig. 7 and Supplementary Tables 8 and 9). Likewise, Gene Ontology (GO) analysis showed enrichment of GO terms related to synapse organization and function in excitatory neurons from APM organoids (Extended Data Fig. 7e and Supplementary Table 10).

To understand whether APM may change neuronal proportions by differentially affecting the proliferation of neural progenitors, we quantified actively cycling cells using Seurat module scores (Methods). APM increased the proportion of cycling aRG cells (binomial linear mixed-effects model (LMEM): n = 7,941 APM cells, 8,808 CDM4 cells; P = 1.2 × 10⁻¹¹, OR = 8.76, 95% CI = 5.25–14.62), whereas CDM4 preferentially enhanced proliferation in interneuron progenitors (binomial LMEM: n = 2,390 APM cells, 8,984 CDM4 cells; P = 3.7 × 10⁻⁵, OR = 0.07, 95% CI = 0.02–0.21), indicating that the two media favour the expansion of different cell populations, even though all cell types were still represented (Extended Data Fig. 6k and Supplementary Table 11).

Analysis of metabolic heath indicated that transcriptional stress signatures were primarily found in progenitor populations, and that the hypoxia signal was largely confined to the organoid core, consistent with previous observations⁶ (Extended Data Fig. 8a–e, Supplementary Table 11 and Supplementary Information). Notably, excitatory neuronal populations under APM conditions showed no increase over time in any stress-related signatures, indicating sustained neuronal health in APM (LMEM; apoptosis: P = 0.24, hypoxia: P = 0.79; glycolysis: P = 0.42; Extended Data Fig. 8f).

To evaluate transcriptional maturation, we again applied the DIALOGUE-based module scoring approach described in Fig. 1d. As in the CDM4 organoids, APM organoids showed strong correlation between time in culture and maturation module score (n = 19 organoids; Pearson’s r = 0.741) (Fig. 3h,i). Importantly, this approach enabled a direct comparison between excitatory neurons in APM- and CDM4-cultured organoids between 4 months and 1 year (the age at which neuronal representation declines in CDM4 datasets). Excitatory neurons in APM organoids scored consistently higher in maturation index (one-sided F test, P = 4.03× 10⁻², linear models) (Fig. 3j); these differences were significant at all timepoints except for 6 months (likelihood ratio test P values, 4 months, 1.56 × 10⁻³; 6 months, 9.77 × 10⁻¹; 9 months, 5.88 × 10⁻²; 12 months, 3.86 × 10⁻⁶; LMEMs) (Fig. 3k).

Finally, we examined whether activity-permissive culture conditions resulted in corresponding changes in the DNA methylome. Comparative WGBS of organoids cultured in CDM4 or APM at matched timepoints (9 months and 1 year; Extended Data Fig. 9a) identified around 900 DMRs between culture conditions (Extended Data Fig. 9b). Notably, these condition-associated DMRs were enriched for neuron-associated and maturation-related genes (Extended Data Fig. 9c,d and Supplementary Table 12). Despite these locus-specific differences, the overall directionality and temporal structure of age-associated methylation remodelling were preserved across culture conditions. Both CDM4- and APM-cultured organoids showed strong correlation between methylation age and time (Extended Data Fig. 9e).

It is unclear whether organoids can build and maintain active neuronal circuits over years in culture. We therefore directly compared neuronal and network activity between culture conditions. We recorded spontaneous extracellular activity using a three-dimensional (3D) multielectrode array (MEA) platform⁹ (Methods) from both APM and CDM4 organoids at 6, 9 and 12 months of differentiation (total of 19 single organoids from three different donor lines; Fig. 4), comparing organoids that were cultured continuously in APM starting from DIV70, to control organoids derived from the same lines that were cultured in conventional medium (CDM4) and transferred to APM 2 weeks before recording as previously described³⁰. For each organoid, we recorded a minimum of 20 min of spontaneous activity and identified spikes post hoc using kilosort³¹, which we optimized for application to organoids (Methods). The majority of recorded organoids under both APM and CDM4 conditions exhibited periodic bursts of highly correlated activity, known as network bursts. Pharmacological treatment verified that the signal corresponded to action potentials and bursting was mediated by synaptic glutamatergic transmission (Methods and Extended Data Fig. 10a,b). Across all three ages, APM organoids displayed a higher spike rate (analysis of variance (ANOVA), P = 2.60 × 10⁻⁴) and increased network burst frequency (ANOVA, P = 4.87 × 10⁻³) (Fig. 4b, Supplementary Table 13, Extended Data Fig. 10c,d and Supplementary Video 3), more spikes per burst (ANOVA, P = 1.25 × 10⁻⁸) and a greater number of neurons participating in each burst (ANOVA, P = 1.76 × 10⁻⁷) (Fig. 4b and Extended Data Fig. 10c,d). Notably, starting from 9 months of differentiation, APM organoids displayed more complex bursting patterns, while CDM4 organoids showed increasingly sparse activity over the same period (Fig. 4c and Extended Data Fig. 10c,d). The most pronounced differences between the two conditions were observed at the 1-year mark: organoids cultured in CDM4 for 1 year consistently failed to exhibit any network bursts (8 out of 8), whereas all APM organoids at 1 year exhibited robust network bursts (9 out of 9).

Notably, recordings from APM organoids that had been maintained in culture for 2 years still showed active bursting, indicating that APM supports the continued presence of active, functional neurons over multiyear timelines (Fig. 4d, Extended Data Fig. 10e and Supplementary Video 4). Notably, the structure of these bursts differed from younger organoids. In organoids that are one year old or younger, when network/population bursts were observed (either in CDM4 or APM), the majority of neurons participate in these events. These neurons fire during every burst, indicating a single, homogeneous mode of collective activation. By contrast, in two-year-old APM-cultured organoids, at least two distinct neuronal subgroups could be identified during bursts, each displaying its own form of collective activation, characterized by unique frequencies, amplitudes and participating neuronal populations (Fig. 4d). In agreement with the functional data, 2-year-old APM-cultured organoids showed the presence of NeuN⁺ neuronal cells by immunohistochemistry (Fig. 4e).

Taken together, these results show that promoting neuronal activity through culture in activity-permissive medium results in the maintenance of functional neuronal populations for at least 2 years in culture. Neurons produced in organoids under these conditions display greater morphological and transcriptional maturity compared with conventional medium, while still being generated according to endogenous molecular programs of development. Notably, our ability to maintain and profile excitatory neurons over these extended culture periods has enabled the generation of a comprehensive multi-modal dataset, comprising molecular, morphological and electrophysiological data on human excitatory cortical neurons at an unprecedented, advanced stage of in vitro development.

Elegant work in the developing mouse and chick embryos has shown that developing cells track developmental time, progressively restricting their fate potential to execute temporally appropriate developmental programs even if transplanted into a heterochronic host^32,33. As a measure of functional changes over time in organoids, we tested whether cells in older organoids retain a memory of the time spent in culture; specifically, whether cells of different ages would show distinctions in fate potential, manifested as their ability to generate early versus late progeny.

To this end, we first tested whether co-development of progenitors of different ages inside the same organoid could be used to probe the differential fate potential of neural progenitors developed for longer or shorter periods of time. For this, we developed an ex vivo mouse chimeroid system using progenitors derived from the embryonic mouse brain (Supplementary Information), as it gives the distinct advantages of fast development and endogenous origin of the progenitor cells. Young organoids (generated from dissociated embryonic day 11.5 (E11.5) mouse cerebral cortex) were cultured for 2 to 3 weeks to generate ‘old’ organoids (Extended Data Fig. 11a–d). We confirmed that neuronal cells were not carried over through dissociation and reaggregation of chimeroids using AAV PHP.eB labelling^34,35,36 (Extended Data Fig. 11e,f), consistent with our previous findings in human cortical chimeroids⁹.

To create chimeroids, we either mixed progenitors from the 2-to-3-week-old organoids alone (monochronic), or mixed cells from the 2-to-3-week-old organoids with freshly dissociated, young, E11.5 cells (heterochronic; Extended Data Fig. 11g). At 1 day after aggregation, old monochronic mouse cultures showed minimal SATB2⁺ cells, while the heterochronic condition exhibited significantly more (ANOVA, P = 6.14 × 10⁻¹⁷, Tukey post hoc P = 2.35 × 10⁻¹⁴) (Extended Data Fig. 11h). Control ‘young’ monochronic chimeroids (derived from E11.5 cortex without previous culture) took longer to begin to produce SATB2⁺ cells, consistent with the fact that they had not yet completed earlier developmental steps (Extended Data Fig. 11i and Supplementary Table 14). Together, these findings indicate that this system can be used to test the fate potential of progenitors.

To test whether human progenitors similarly change their fate potential as a reflection of the time passed in culture, we produced human chimeroids⁹ by mixing either old cells (from organoids collected at 9–12 months), or a combination of old and young (from organoids at DIV15) cells.

First, we produced old monochronic chimeroids (monochronic (old)) from late-stage cortical organoids (9–12 months old) (Fig. 5a,b), and profiled pooled chimeroids (15–20 chimeroids per batch) at DIV15–19 after reaggregation using scRNA-seq. The monochronic (old) samples mirrored the cellular composition expected for organoids grown for 9 months of culture, containing all major cell types (over 5% abundance) found in standard organoids at 9 months. Of these major cell types, the monochronic (old) samples had a reduced prevalence of glial precursors (negative binomial generalized LMEM, P = 8.5 × 10⁻⁴) and oRGs (P = 0.038); the abundances of immature INs, IN progenitors, tRGs, aRGs and astrocytes were not significantly different (P > 0.1), despite the monochronic (old) chimeroids having developed for only 15 days after reaggregation (Fig. 5a,c,d and Supplementary Tables 5 and 14). Both were distinct from the cell types made in young chimeroids at DIV15 (aRGs, IPs, subplate cells, cortical hem cells, subcortical populations, Cajal Retzius cells and choroid plexus cells)⁹. This suggests that old progenitors retain a memory of the time they had previously spent in culture, producing progeny that is appropriate for their age (9–12 months), such as astroglia.

It is possible that the continued production of late cell types in the monochronic (old) chimeroids may be the result of a lack of signals instructive of earlier fates within these chimeroids. We therefore tested the fate potential of these old progenitors when exposed to instructive conditions from early cell fates by co-culture with cells from younger organoids (heterochronic chimeroids, or heterochronic). Progenitors from late-stage cortical organoids (9 to 12 months old) were mixed with developmentally asynchronous progenitors from young organoids (DIV15), and profiled using scRNA-seq at DIV15 after aggregation (Fig. 5f–h). To resolve the originating organoid (and therefore age) contributing to each differentiated cell type in the chimeroids, we used different donor cell lines for young (PGP1) and old (GM) organoids^9,37 (Methods). Notably, old progenitors within heterochronic cultures (that is, those derived from the 9 month organoids) exhibited lower maturation scores than progenitors from monochronic (old) cultures (Extended Data Fig. 11j). Similarly, GSEA of genes that change during ageing from 6 to 9 months in canonical cortical organoids showed that old progenitors in heterochronic organoids shared a greater proportion of molecular programs with 6-month progenitors than monochronic (old) progenitors did (Extended Data Fig. 11k). Together, these results indicate that the younger cells in the heterochronic cultures provided instructive signals that the older cells are able to respond to.

Notably, derivatives of the old cells in heterochronic organoids showed a substantially higher proportion of CPNs (49.0%) compared with both the monochronic (old) group (replicate 1 = 0.55%, replicate 2 = 0.50%, replicate 3 (11 months) = 0.02%) and 9-month standard organoids (CPNs; 1.1%) (Fig. 5h,i). In agreement, while both old (9 month) and young (1 month) monochronic cultures had few or no putative CPNs (SATB2⁺ cells), SATB2⁺ putative CPNs were clearly detected in the heterochronic condition by immunohistochemistry (all assayed at DIV15 after reaggregation; Fig. 5j). These results indicate that exposure to inductive signals from young progenitors induced the old progenitors to restart excitatory neurogenesis, which occurs before astrogliogenesis. Notably, while the young-derived cells within the same heterochronic chimeroids produced only the early cell types expected at DIV30 (15 DIV + 15 DIV after aggregation), including cortical hem cells, subcortical precursors, aRGs, IPs, newborn deep-layer neurons and choroid plexus, the old progenitors produced CPNs, glial precursors and astrocytes, late-derived populations that are normally produced at 2–3 months of culture. Thus, in only 2 weeks after reaggregation, the old progenitors can generate late-stage neurons (CPNs), while their young progenitor counterparts produce much earlier cell fates within the same chimeric organoid.

Together, the human and mouse results indicate that cortical organoid cells retain a form of temporal memory. In human organoids, old progenitors are still plastic after 9 months in culture and can respond to instructive neurogenic signals to make excitatory neurons again; however, they have recorded the passage of time and retain a memory of the developmental steps already performed, such that even when exposed to signals instructive of an earlier neurogenic fate, they ‘skip’ production of early progeny to generate excitatory neurons normally produced after months in culture.

The speed of development of cells and organs is markedly different across species³. The human brain displays extremely long periods of neoteny, taking close to two decades to form and mature. Human brain organoids that can continue to develop over years, rather than weeks and months, offer a unique opportunity to experimentally access species-specific processes of human brain postnatal maturation that cannot be easily accessed in vivo.

By enabling culture of human brain organoids for over 5 years, we demonstrate that organoids develop through similar transcriptional and epigenetic programs associated with age progression in the human brain in vivo, eventually acquiring features of human postnatal ages. These findings demonstrate that temporally structured regulation of cortical development can unfold over unprecedented timescales, even outside the context of the body.

It is notable that, during this process, organoids record the passage of time, and can also recall the time spent in culture, as shown by the fact that cells derived from old organoids can skip steps of development that they have already performed, and generate in 2 weeks late neuronal progeny that would normally take over 2 months to produce.

The data provide additional motivation for efforts to further integrate into organoid systems additional driving forces that shape brain development in vivo. Our demonstration that promoting spontaneous activity in organoids can support neuronal health and maturation over multiple years points at the potential benefit of incorporating factors that have been demonstrated to be critical for brain development in vivo, such as evoked activity through sensory systems and the influence of non-neural tissues.

Together, organoids cultured over extended timelines and the multimodal wealth of data produced from them represent both a powerful experimental system and the information to fuel understanding of the largely unexplored mechanisms governing human brain neoteny, maturation and evolution.

Ethics statement

All experiments involving human cell lines were approved by the Harvard University IRB and ESCRO committees. All animal experiments were conducted according to protocols approved by the Institutional Animal Care and Use Committee (IACUC) of Harvard University. All experiments were performed in accordance with relevant guidelines and regulations, and in accordance with informed consent obtained from the donors of the originating cells or tissues.

Mouse experiments and housing

To develop the mouse-endogenous chimeroid system, we used wild-type CD-1 and C57BL/6 mice (Charles River Laboratories). We used embryos of both sexes, collected at E11.5. The sex of each embryo was not determined prior to the experiment. Sample size was not predetermined using specific methods. Mice were maintained under standard housing conditions in a temperature- and humidity-controlled facility under a 12 h–12 h light–dark cycle (light, 07:00 to 19:00) with ad libitum access to food and water. The ambient temperature was maintained at 20–24 °C and the relative humidity at 30–70%.

Human pluripotent stem cell culture

All human PS cell lines were maintained as previously detailed^7,9. In brief, MTESR1 medium (StemCell Technologies), mTESR⁺ medium (StemCell Technologies) or StemFlex medium (Gibco), all with 1% of added penicillin–streptomycin Solution (Corning), were employed for culture of stem cells in cell culture dishes (Falcon) precoated with 1% Geltrex (Gibco), at 37 °C in 5% CO₂. All human PS cells were maintained below passage 55 and tested negative for mycoplasma (assayed using the MycoAlert PLUS Mycoplasma Detection Kit, Lonza).

Characterization of the PS cell lines

The psychiatric control Mito210 male iPS cell line was provided by B. Cohen (McLean Hospital); the PGP1 male iPS cell line was provided by G. Church; the GM08330 male iPS cell line was provided by M. Talkowski (MGH) and was originally derived from fibroblasts obtained from the Coriell Institute for Medical Research^6,7,30. The H1 male human embryonic stem cell line (also known as WA01) was purchased from WiCell; and the 11a iPS cell line was obtained from the Harvard Stem Cell Institute. The female CW50037 iPS cell line was from the California Institute for Regenerative Medicine (CIRM) iPS cell collection. All lines were authenticated as follows. The PGP1 iPS cell line was authenticated by short-tandem-repeat (STR) analysis (performed by TRIPath). The 11a cell line was authenticated by karyotyping as previously described⁴⁰. The Mito210 iPS cell line was authenticated with genotyping analysis (Fluidigm FPV5 chip) performed by the Broad Institute Genomics Platform. The CW50037 iPS cell line was authenticated using single-nucleotide polymorphism (SNP) genotyping (Illumina Global Screening Array (GSA) from Illumina; processed at the Genomics Platform, Broad Institute) for cell line identification and detection of chromosomal abnormalities. The H1 and GM08330 lines were authenticated by STR analysis (performed by WiCell). The GM08330 parental line has a previously reported³⁰ interstitial duplication in the long (q) arm of chromosome 20; all other lines were karyotypically normal. No commonly misidentified lines were used in this study.

Cortical organoid and chimeroid differentiation

Dorsally patterned cortical organoids were generated following a protocol previously published by our group⁷. In brief, on day 0, feeder-free cultured human PS cells, 75–85% confluent, were enzymatically dissociated to single cells with Accutase (Gibco), and 9,000 cells per well were reaggregated in ultra-low-cell-adhesion 96-well plates with V-bottom conical wells (sBio PrimeSurface plate; Sumitomo Bakelite) using the same pluripotent cell medium in which they were previously maintained. At day 1, 80 μl of medium was replaced with cortical differentiation medium (CDM) I, containing Glasgow-MEM (Gibco), 20% knockout serum replacement (Gibco), 0.1 mM minimum essential medium non-essential amino acids (MEM-NEAA) (Gibco), 1 mM pyruvate (Gibco), 0.1 mM 2-mercaptoethanol (Gibco), with 1% of penicillin–streptomycin solution (Corning). From day 0 to day 6, ROCK inhibitor Y-27632 (Millipore) was added to the medium at a final concentration of 20 μM. Patterning small molecules, WNT inhibitor IWR1 (Calbiochem) and TGFβ inhibitor SB431542 (Stem Cell Technologies), were added from day 0 to day 18, at a concentration of 3 μM and 5 μM, respectively.

Single-donor and multidonor chimeroids were generated as previously described⁹. In brief, patterned EBs between day 15–18 were dissociated into single-cell suspensions using a modified papain-based protocol⁹. After confirming morphology, organoids were enzymatically and mechanically dissociated, then filtered and resuspended in CDM1 medium with ROCK inhibitor. Cell suspensions from different donors (in the case of multidonor chimeroids) were mixed in equal ratios and reaggregated (18,000–20,000 cells per well) in ultra-low-adhesion 96-well plates. Then, 2 days later, embryoid bodies were transferred to low-attachment dishes and cultured under orbital agitation. Medium was sequentially changed to CDM2, CDM3 (at DIV35) and CDM4 (at DIV70) as per established protocols⁷. For APM, organoids at day 70 were transferred into a medium containing BrainPhys Basal (Stem Cell Technologies), 20 ng ml⁻¹ N2 supplement (1×, Thermo Fisher Scientific), 20 ng ml⁻¹ B27 supplement (1×, Thermo Fisher Scientific), 1% penicillin–streptomycin solution 100× (Corning), 1% MEM-NEAA (Gibco), 1% GlutaMax (Gibco) and amphotericin B (Thermo Fischer Scientific). After filtration, GDNF and BDNF at a final concentration of 20 ng ml⁻¹ were added. Along with dCAMP (MilliporeSigma), ascorbic acid (StemCell Technologies) and laminin (Thermo Fischer Scientific) at a concentration of 1 mM, 200 nM and 1 μg ml⁻¹. Half of the medium was replaced with fresh medium twice per week. Heterochronic and monochronic cultures were obtained using the previously described chimeroid strategy; however, papain dissociation time and mechanical dissociation were adjusted based on organoid age (25 min for younger organoids; for older organoids, mincing with a blade followed by 45 min of digestion). While the chimeroid strategy is extremely efficient in younger organoids, the success rate of the heterochronic cultures remains very low. This is probably due to the dual requirement for sufficient patterning in the young organoids and adequate stemness in the older cells, both of which are necessary to ensure robust signalling from the young tissue and sufficient recovery of older cells for analysis at day 14 after mixing. Chimeroids derived from old organoids were smaller than those derived from young organoids, probably due to the limited proliferative potential of the old neural progenitors. Organoid time-series include at least n = 6 organoids/3 batches at all stages (6 months to 5 years).

Mouse experiments

We used wild-type CD1 and C57BL/6 mice (Charles River Laboratories). Pregnant dams were euthanized to obtain embryos. Embryos were staged carefully according to Theiler stages for E11.5 and the neocortex was dissected. Stage-matched embryos were pooled (15–20 per experiment) prior to dissociation. The dissociation protocol described for the human chimeroids was used with some modifications. Dissociated progenitors were seeded to aggregate overnight at a density of 12,000 cells per well. The next day, the medium was switched to CDM2 for 3 days. On day 4 after seeding, the medium was changed to CDM3 and then to CDM4 on day 6. We quantified the number of nuclei with positive SATB2 expression across three condition groups: monochronic (old), monochronic (young) and heterochronic from cultures 1 day (D1AM; for each group: n = 3 replicate organoids × 3 batches) and 2.5 days (D2.5AM; for each group: n = 3 replicate organoids × 1 batch) after mixing, respectively. For evaluating the overall effect of condition on the proportion of nuclei with positive SABT2 expression at D1AM, binomial generalized linear models (GLMs) were used to model the number of SATB2-positive nuclei out of total nuclei per organoid and adjust for batch by including batch as a fixed effect. Models with or without adding a condition as a fixed effect were evaluated by a likelihood ratio test using the anova function. Following the observation of the global effect of the condition being statistically significant (ANOVA, P < 0.05), estimated marginal means for each group were derived from the fitted full binomial GLM by using the R package emmeans, and Tukey-adjusted pairwise comparisons were performed between groups. For D2.5AM, due to the absence of more than one batch level to adjust for, we proceeded with computing estimated marginal means and performing pairwise comparisons between groups.

To fluorescently label neurons, we incubated cells while aggregating with AAV particles encoding GFP under the CAG promoter. The AAV serotype PHPeB permits efficient transduction of the central nervous system³⁶. Organoids were fixed and imaged around 3weeks after aggregation on the LSM900 microscope. Images were acquired using the Zeiss microscope and processed in ImageJ (v.2.14.0/1.54 f). pAAV-CAG-GFP was acquired from Addgene (Addgene, 37825).

Statistical analysis

Organoids used for analysis or treatment were chosen from each batch to be representative of the morphology seen in that differentiation. For every experiment, organoids were collected without a preconceived selection strategy or priority. The investigators did not use blinding in this study. However, our analytical pipeline for each experiment followed uniform criteria applied to all samples, allowing us to analyse our data in an unbiased manner. All bioinformatics analyses were applied uniformly to all samples without considering genotype adjustments. Sample size was not predetermined using specific methods. At least three organoids were used for each experiment, guided by analysis from previous published work^6,7,30. This prior research demonstrated the reproducibility of this organoid system and confirmed that using three organoids adequately captured variation.

Fixation and processing of samples for cryosectioning

Organoids were fixed in 4% paraformaldehyde (PFA) (Electron Microscopy Services) overnight in a 12-well plate (Falcon) at 4 °C, washed three times with 1× PBS (Gibco), and cryoprotected in a 30% sucrose solution (Sigma-Aldrich) in PBS overnight at 4 °C.

Gelatin solution containing 10% bovine gelatin (Sigma-Aldrich) and 7.5% sucrose (Sigma-Aldrich) was prewarmed at 37 °C for 15 min. The 30% sucrose solution was removed from the samples and exchanged for the prewarmed gelatin, and the samples were incubated at 37 °C for 15 min. Meanwhile, plastic moulds were coated with a 2 mm layer of warm gelatin solution and left to polymerize at room temperature. The samples were then transferred to the pretreated plastic moulds, and 1 ml of warm gelatin solution was added on top. After polymerization at room temperature for 3 min, the samples were prechilled at 4 °C for 15–20 min. Finally, the moulds were frozen in a cold bath containing 100% ethanol and dry ice for 2–3 min, and stored at −80 °C indefinitely.

Immunohistochemistry

For immunohistochemistry, 14–18 μm-thick sections were cut using the cryostat (Leica). Cryosections were stabilized at room temperature for 5 min and blocked with 10% donkey serum (Sigma-Aldrich) + 0.3% Triton X-100 (Sigma-Aldrich) in PBS. Primary antibodies (Supplementary Table 15) were diluted in the blocking solution and incubated overnight. After four washes with PBS, cryosections were incubated at room temperature with secondary antibodies diluted in PBS (1:1,000; Supplementary Table 15) for 1 h at room temperature, washed four times with PBS and stained with DAPI (1:10,000 in PBS + 0.1% Tween-20) for 5 min to visualize cell nuclei. At later timepoints (>1 year), we consistently observed a reduction in staining specificity, with increased background signal that limited reliable interpretation of most antibodies, despite extensive optimization. For the 2–3 year organoid immunohistochemistry, we used the fresh-frozen protocol described below.

Immunohistochemistry on fresh-frozen sections

Human brain organoids aged 2 to 5.8 years were rinsed with PBS, directly embedded in OCT compound (NEG-50, Richard-Allen Scientific) and then quickly frozen. Fresh-frozen OCT blocks were cryosectioned at 10 µm thickness. The sections were thawed onto Superfrost Plus Gold glass slides (EMS). Organoid depths of 0–150 μm were collected as serial sections. Sections were fixed in ice-cold methanol at −20 °C for 20 min, then rinsed with PBST0.1 (1× PBS with 0.1% Tween-20). Permeabilization was performed with PBST0.25 (1× PBS with 0.25% Tween-20) for 15 min. Blocking was done with 5% donkey or goat serum in PBST0.5 (1× PBS with 0.5% Tween-20) for 30 min. Alexa-Fluor-conjugated antibodies (NEUN/RBFOX3; 608455, BioLegend) were diluted 1:200 in blocking buffer and incubated with sections for 1–2 h. The sections were rinsed with PBST0.1 and incubated with 0.2 μg ml⁻¹ DAPI for 3 min, then rinsed with PBS and water. ProLong Gold anti-fade mounting medium was applied, and coverslips were mounted. Images were acquired using LSM900 with Zeiss software. Image processing was done using ImageJ (v.2.14.0/1.54f). NeuN-positive nuclei were observed up to approximately 80 μm deep, beyond which they were not observed.

Hypoxyprobe assay

Hypoxyprobe (pimonidazole HCl; Hypoxyprobe Kit HPI-100) was applied to n = 3 1-year-old CDM4 organoids at a final concentration of 100 µM. Organoids were incubated with the probe for 1.5 h at 37 °C and 5% CO₂. Samples were then washed with PBS, fixed in 4% PFA at room temperature for 2 h, embedded and sectioned at 14 µm. Immunohistochemistry was performed according to the manufacturer’s instructions, and the primary antibody included in the kit was used at a 1:50 dilution⁶.

Whole-organoid immunofluorescence

Organoids were washed in wash buffer (PBS with 0.4% BSA) and fixed at room temperature for 30 min with 4% PFA. After fixation, organoids were stored in 1× PBS with 0.1% Tween-20 (P9416, Sigma-Aldrich). Permeabilization, blocking and staining were performed as previously described⁴¹. Nuclei were counterstained with 0.5 μg ml⁻¹ DAPI for 30 min. Organoids were optically cleared using RIMS and mounted as previously described^41,42. Images were acquired with LSM880 (Zeiss) and LSM900 (Zeiss) using a ×40 objective. Images were processed with ImageJ (v.2.14.0/1.54f). Nuclei were segmented and quantified using ZEN blue v.3.1 Image Analysis and Intellesis software packages. Image segmentation was achieved using a previously trained ZEN Intellesis algorithm⁴³. SATB2 intensity was thresholded to define positive nuclei and counted across the various z-stacks of individual organoids.

Microscopy

Immunofluorescence images were acquired using the Zeiss Axio Imager.Z2 with a ×20 objective (pixel size, 0.325 μm) using the Apotome optical sectioning function. The Zen Blue software was used to perform tile stitching and apotome deconvolution before exporting the images. For the LSM900, z stacks with 3 μm step size were acquired with a ×20 objective (pixel size, 0.62 μm), followed by tile stitching using the Zen Blue software. Further processing, such as z projections, channel merging and adding scale bars, was performed in Fiji v.43. Images were adjusted for brightness and contrast; all adjustments were applied to the whole image.

SATB2 and FOS quantification

Nine-month-old organoids from three genetic backgrounds (H1, PGP1 and 11a) were analysed (n = 3 organoids per condition per line; APM and CDM4). Three regions per organoid were randomly selected based on DAPI signal and imaged from three sections on the same slide. z-Stacks were acquired using a 20× objective (0.62 µm per pixel; 3 µm step size; 3–5 optical sections) on a Zeiss system (Zen Blue). Acquisition settings were as follows: SATB2 (650 V, 2%, 28 µm pinhole), FOS (750 V, 10%, 31 µm) and DAPI (650 V, 1%, 27 µm). Stacks were converted to average-intensity projections and quantified in Fiji (v.43). SATB2-positive cells were defined as DAPI-positive nuclei with elevated SATB2 signal and were subsequently assessed for FOS expression.

Bulk RNA-seq from organoids

Snap-frozen organoids were resuspended in RLT buffer containing β-mercaptoethanol, and RNA was isolated with a DNase digestion step according to the manufacturer’s protocol (RNeasy, Qiagen). Then, 10 ng of RNA was used for library preparation with the Smart-Seq Pico input total RNA library preparation kit, including ZapR depletion (with unique molecular identifiers (UMIs)), according to the recommended protocol. Libraries were quantified, pooled and sequenced on the NextSeq 2000 (Illumina) instrument.

DNA extraction and methylome analyses

Genomic DNA extraction and WGBS

Genomic DNA was isolated from organoids using a previously described method⁴⁴. In brief, organoids were thawed on ice and suspended in 500 μl of resuspension buffer (240 μl of genomic lysis buffer (D4075, Zymo), 20 μl of proteinase K (D4075, Zymo) and 240 μl of nuclease-free water). The sample was thoroughly mixed, vortexed and incubated at 55 °C for 4 to 7 h at 300 rpm. An equal volume of phenol:chloroform:isoamyl alcohol (15593031, Thermo Fisher Scientific) was added, mixed by inversion and then centrifuged at 13,000g for 5 min. The aqueous phase (around 400 μl) was transferred to a new tube for precipitation. To this aqueous phase, 8 μl of Glycogen-Blue (AM9515, Thermo Fisher Scientific), 20 μl of 5 M NaCl (AM9760G, Thermo Fisher Scientific) and 1 ml of ethanol (E7023, Sigma-Aldrich) were added, mixed by inversion, and allowed to precipitate overnight at −20 °C. The tubes were spun at 13,000g for 45 min at 4 °C. The supernatant was discarded and 1 ml of 70% ethanol was added to the pellet, which was then spun down at 13,000g for 45 min at 4 °C. The supernatant was discarded, and the pellet was air-dried for 10 min before being resuspended in elution buffer (low-EDTA TE pH 8, 15575020, Thermo Fisher Scientific) and incubated at 55 °C for 10 min.

Genomic DNA was fragmented in a tube (MicroTube AFA Fiber pre-slit snap-cap 6 × 16 mm, 520045, Covaris) using a Covaris S-series S2 Sonicator with the following settings: duty cycle, 10%; intensity, 5; cycles per burst, 200; maximum temperature, 7 °C; duration, 76–90 s. Fragmented DNA was purified and concentrated using the DNA Clean & Concentrator kit (D4013, Zymo) and eluted in 22 µl of low TE (pH 8). 1 µl of the eluate was used for quality control on the Agilent TapeStation (HSD5000), resulting in an average fragment size of 237–322 bp. Bisulfite conversion was performed using the EZ DNA Methylation-Gold kit (D5005, Zymo), with the final converted DNA eluted in 16 µl of low TE (pH 8). Libraries were prepared using the xGen-Methyl-Seq DNA Library prep kit (10009824, IDT), in combination with xGen UDI Primer Plate 2 (8 nucleotides, 10009816, IDT). The number of PCR cycles used was 7. Each library underwent two final rounds of purification using AMPure XP beads (A63881, Beckman Coulter). Libraries were assessed for concentration, fragment size distribution and primer-dimer absence using the Agilent TapeStation (HSD5000). Sequencing was conducted on the Illumina NovaSeq 6000 and Element Biosciences Aviti platform, generating 150 base paired-end reads, targeting 800–900 million reads per sample.

WGBS processing

Raw reads were subjected to adapter and quality trimming using cutadapt⁴⁵ (v.4.6; parameters: --quality-cutoff 20 --overlap 5 --minimum-length 25; Illumina TruSeq adapter clipped from both reads), followed by trimming of 10 and 5 nucleotides from the 5′ and 3′ end of the first read and 15 and 5 nucleotides from the 5′ and 3′ end of the second read. Trimmed reads were aligned to the human reference genome (hg19) using BSMAP⁴⁶ (v.2.90; parameters: -v 0.1 -s 16 -q 20 -w 100 -S 1 -u -R). In this study, we used hg19 (rather than the current hg38) as the reference genome to maintain compatibility and coordinate consistency with previously generated reference datasets, annotation resources and downstream comparative analyses incorporated into our methylation analysis pipeline, including publicly available developmental brain methylation datasets and established epigenetic clock resources. Importantly, all samples and comparative analyses were processed uniformly within the same reference framework. Because our analyses focused primarily on genome-wide and regional DNA methylation patterns rather than fine-resolution structural variation or novel genomic annotations, we do not expect the use of hg19 versus hg38 to materially affect the biological conclusions of the study. Sorted BAM files were generated and indexed using samtools⁴⁷ with the ‘sort’ and ‘index’ commands (v.1.21). Duplicates were removed using the MarkDuplicates command from GATK (v.4.5.0.0; --VALIDATION_STRINGENCY = LENIENT --REMOVE_DUPLICATES=true --COMPRESSION_LEVEL 4 --ASSUME_SORT_ORDER coordinate)⁴⁸. Methylation rates were called using mcall from the MOABS⁴⁹ package (v.1.3.9.6; default parameters; --reportCpX A/C/T for non-CpG methylation calling). All methylation analyses were restricted to autosomes and only CpGs covered by at least 10 and at most 150 reads were considered for downstream analyses. Assessment of global, genome-wide methylation was performed using average (arithmetic mean) methylation levels. Region-specific methylation was assessed by calculating average methylation levels across features (bins, partially methylated domains, highly methylated domains, DMVs, repeats, cdDMRs, superenhancers) using bigWigAverageOverBed from UCSC tools for each sample, where a feature was only considered if at least three CpGs were contained measurements within a region.

Endogenous comparison

cdDMRs defined from postnatal human brain were obtained from ref. ¹⁴. cdDMR coordinates in hg19 and their cluster assignment based on cell-type-specific methylation trajectories (six clusters: 1:G−N+, 2:G0N+, 3:G0N−, 4:G+N0, 5:G+N−, 6:G−N0) were used as provided by the study. cdDMRs were scored for mean methylation in organoid samples using bigWigAverageOverBed as described above (≥3 CpGs covered). To quantify recapitulation of temporal methylation dynamics in organoids, for each cdDMR, a linear regression model was fitted predicting mean cdDMR methylation from time in culture (in months), yielding a slope estimate (change in methylation per month). cdDMRs were classified as showing temporal changes over time in organoids if the absolute slope value had a minimum of 0.1/60 (corresponding to a predicted minimum change of 0.1 in methylation rate over 60 months). The fetal brain DNAm age was estimated using the FetalClock function from ref. ¹⁸, based on the implementation provided at GitHub (https://github.com/LSteg/EpigeneticFetalClock).

Differential methylation analysis

DMRs were called using metilene⁵⁰ (v.0.2-8). DMRs were defined by an absolute minimum difference in methylation of 0.1 with a maximum distance of 300 nucleotides between CpGs within a DMR and a minimum of 10 CpGs per DMR and Bonferroni correction for multiple testing (parameters: --m 10 --c 1 --d 0.1) and filtered by q < 0.05. DMRs were classified as hypo-DMRs and hyper-DMRs based on the direction of change. DMRs were associated with CGIs and DMVs based on a minimum overlap of 1 bp. DMRs were annotated to the nearest genes using GREAT⁵¹ (v.4.0.4; default parameters).

To identify temporal changes in methylation, DMRs were called between early and late culture timepoints (3 months versus 5 years). DMRs were filtered to retain only those with a minimum absolute Pearson correlation coefficient of 0.6 between time in culture (in months) and mean DMR methylation rate across all timepoints.

To evaluate effects of APM on the methylation landscape, DMRs were called between organoids cultured for 9 months or 1 year in CDM4 and organoids cultured for 9 months or 1 year in APM.

Genomic feature annotation

CGIs were obtained from the UCSC Genome Browser (https://genome-euro.ucsc.edu/cgi-bin/hgTables). CpG shores were defined as 2 kb regions flanking CGIs upstream and downstream; shelves were defined as 2 kb regions flanking the shores.

DMVs were defined based on WGBS data from the HUES64 human embryonic stem cell line⁵². To identify DMVs, a previously described sliding-window approach was adapted⁵³. In brief, 5 kb windows with a 1 kb step size were generated across autosomes using bedtools⁵⁴ makewindows. Average CpG methylation was calculated for each window, excluding CpGs located within CGIs, and separately for CGIs. Windows and CGIs with an average methylation below 0.15 and containing at least 10 CpGs were merged, excluding regions composed solely of CGIs.

Annotations of highly methylated domains and partially methylated domains were obtained from GitHub (https://zwdzwd.github.io/pmd)⁵⁵.

Repeat annotations were obtained from the hg19 UCSC RepeatMasker track, excluding entries on alternative haplotypes, fix patches and chrMT.

Epigenetic clocks

CpG methylation values at low-coverage sites (coverage < 5) were imputed using the boostme package (v.0.1.0). Imputation was performed on individual samples represented as objects from the bsseq package (v.1.40.0). Imputed values exceeding one were capped at one. Array-based probe IDs were mapped to genomic coordinates using the Illumina EPIC array manifest (EPIC.hg19.manifest). The pan-tissue human methylation clock¹⁷ was applied using the DNAmAge function from the methylclock⁵⁶ package (v.1.10.0). Cortical DNA methylation age⁵⁷ was estimated using the CorticalClock function based on the implementation provided at GitHub (https://github.com/gemmashireby/CorticalClock). DNAm age was assessed by the median absolute error between the time in culture of the organoids and the predicted methylation age in months, as well as Pearson correlation between the two across all samples.

To assess mitotic history, DNA methylation at solo-WCGW CpGs, which are particularly susceptible to methylation loss with cell division⁵⁵, was quantified. Methylation was evaluated across three feature sets: all solo-WCGWs, solo-WCGWs within common PMDs and solo-WCGWs in common PMDs represented on the HM450 array. Annotations were obtained from https://zwdzwd.github.io/pmd, and the average methylation levels across feature sets were calculated per sample from non-imputed WGBS data.

mCA

Forebrain superenhancers were obtained from a previously published dataset of fetal human brain tissue and cerebral organoids¹⁵. The mean mCA levels at these regions were compared with the background levels calculated across genome-wide 30 kb bins (approximating the average superenhancer size), excluding bins overlapping with superenhancers.

Locus methylation trend plots

For regions of interest (ROIs), the mean methylation was computed using bigWigAverageOverBed as described above (≥3 CpGs covered) and modelled as a function of time in culture (in months) using linear regression. Trends were visualized using per-ROI scatterplots with fitted regression lines (and 95% confidence intervals) and Pearson correlation estimates.

Visualizations

Unless stated otherwise, all statistics and plots were generated using R v.4.4.1. Violin plots were generated using the vioplot (v.0.5.1) or ggplot2 (v.3.5.2) package and show the kernel density estimation with embedded box plots indicating the median, interquartile range and whiskers extending to 1.5× the interquartile range. Circos plots were created using the RCircos (v.1.2.2) package from tracks with averaged profiles over 50 kb bins created using UCSC tools bigwigAverageOverBed. The differential signal was obtained by subtracting the 5-year profile from the 3-month profile. Cytoband and ideogram data for hg19 were obtained from the RCircos package. DMR methylation heat maps and average tracks were created using the EnrichedHeatmap⁵⁸ package (v.1.34.0) by applying the normalizeToMatrix function first, with 2,500 bp extensions up- and downstream and a target ratio of 0.4. Pairwise genome-wide correlations of methylation rates between samples of consecutive timepoints were plotted using the smoothScatter function. The box plots show the median (centre line), interquartile range (box) and the whiskers extend to 1.5× the interquartile range; individual observations are overlaid as points. Bar plots, box plots, scatter plots and line plots were created using ggplot2.

Electrophysiological analyses

Electrical signals were recorded using the Accura 3D CMOS-HD-MEA system by 3Brain (BioCAM DupleX system, in combination with Accura 3D chips). The 3D chip contains 4,096 penetrating µNeedle electrodes with 90 µm height arranged in 64 × 64 grid in a square measuring 3.8 × 3.8 mm, with a sampling rate of 20 kHz and a 12-bit resolution. CDM4 organoids were changed to APM medium 14–21 days before recording, and acute recordings from intact organoids were performed at 37 °C in a mini-incubator box filled with Carbogen. BrainWave v.5 was used to record spontaneous activity for 15–20 min. To test whether the electrical activity measured was synaptic, NMDA and AMPA receptor activity was blocked at the end of recordings using D-AP5 (150 μM) and DNQX (30 μM), respectively. Moreover, action potentials were blocked by bath application of TTX (1 μM). Spike detection was carried out by Kilosort2³¹, followed by analysis of the spike train using custom MATLAB scripts (v.24.1, R2024b, The MathWorks). The detection of rapid spiking periods (burst) was performed by implementing the maximum interval⁵⁹ algorithm (maximum interval = 170 ms, maximum end interval = 400 ms, minimum interval between bursts = 800 ms, minimum duration of burst = 40 ms, minimum number of spikes = 4). The analysis of network bursting was performed on the basis of the population-averaged spiking rate along all detected units. A peak in the population signal was considered to be a network burst if it met the following criteria: (1) the peak amplitude was greater than 4× the s.d. of the noise value; (2) a set of bursting cells composed of at least 20% of total cells were active during that population spike; and (3) a cell was considered part of the set of bursting cells only if it participated in at least 50% of the network bursts. The peaks of the network bursts were used to measure the interburst interval (IBI), and the network burst rate was obtained from the average IBI. For organoids that displayed no bursts or networks bursts, the amplitude and duration values were imputed as 0. Outlier datapoints were identified and removed using the ROUT method (GraphPad Prism 10.5.0) with Q = 0.1%.

EM analyses

A total of 18 organoids (three organoids per condition at 6 months: two H1 and one 11a genetic background, six organoids per condition at 12 months for both CDM4- and APM-treated conditions) were immersion-fixed in 2% PFA (EMS, 15710), 2.5% glutaraldehyde (EMS, 16220) and 0.003 M CaCl₂ (Sigma-Aldrich, C4901) in 0.1 M sodium cacodylate buffer (caco buffer, Sigma-Aldrich, C0250) for 48 h (at room temperature for the first 24 h and at 4 C for the remaining 24 h). Organoids were prepared for EM by treating them with 1% OsO₄ (EMS, 19170) with 10 mg ml⁻¹ potassium ferrocyanide (EMS, 20150) and 0.003 M CaCl₂ in 0.1 M caco buffer for 1 h at room temperature. Organoids were washed afterwards and stained with a 2% OsO₄ solution with 0.003 M CaCl₂ in 0.1 M caco buffer for 1 h at room temperature. The sections were then stained with 2% uranyl acetate (EMS, 8473) at 4 °C overnight, dehydrated and embedded in LX-112 epoxy resin (LADD Research Industries, 21310) at 60 °C for 72 h. Blocks were trimmed and sectioned at 40 nm slice thickness with a Leica EM UC6 Ultramicrotome and sections were collected on carbon-coated Kapton tape using an automatic tape-collecting ultramicrotome⁶⁰. Strips of tape were mounted onto silicon wafers and sections were post-stained with 2% uranyl acetate in water and 3% lead citrate (Leica Biosystems ULTROSTAIN II). Sections were imaged using a FEI Magellan scanning electron microscope (Thermo Fisher Scientific) equipped with a custom image acquisition software (WaferMapper)⁶¹. A total of 12 panoramic high-resolution images (one per organoid) were acquired using the backscattered electron detector (7 kV, 1 μs per pixel dwell time, ranging from 86,000 to 106,000 pixels in the x axis and from 106,000 to 150,000 pixels in the y axis, at 4 nm resolution). Owing to natural size variation among organoids, the final 2D panoramas varied slightly in dimensions to ensure spatial coverage from the surface to the central region of each organoid. Stitched EM images were imported into VAST⁶² for visualization and manual annotation of synaptic clefts, creating a ground-truth dataset. A deep neural network based on a U-Net architecture⁶³ was then trained on this ground truth using the mEMbrain MATLAB package²¹. The annotated dataset comprised less than 5% of the total volume and was produced by two human experts (around 20 h total annotation time; pixel-level validation accuracy > 0.95). mEMbrain’s interactive environment facilitated both model training and accuracy evaluation. Predicted synapses were represented as 2D probability maps registered to the EM images in VAST. Synapses were initially defined as connected components of pixels with a synapse probability greater than 40%, a threshold chosen to approximate expert annotations. Model performance exceeded 95% accuracy based on manual verification of the largest synapses (>95th percentile by area), providing an estimate of the false-positive rate for unbiased condition comparisons. To standardize automatic synapse quantification across sections, a computational synapse-size threshold was determined for each section individually. Specifically, 103 randomly selected predicted synapses per section were manually classified as synapse or non-synapse by a human expert. A psychometric curve was then fitted to these annotations to establish a threshold corresponding to over 80% model accuracy. Final synapse counts per section were obtained using these individualized thresholds, ensuring objectivity and consistency across samples. To account for structural variability and region-specific differences in synapse distribution, each section was divided into five horizontal segments. Synaptic density was calculated in each segment, and the one with the highest density was selected for statistical comparison (within each section, synaptic density variation across segments ranged from 0.00 to 0.14). This approach minimized sampling bias across organoid layers. The same set of 103 randomly selected predictions was further used to classify synaptic contacts by their postsynaptic target. Synapses were categorized as spine synapses if they were located on putative dendritic spines—morphologically consistent with spine-like structures observed in 2D EM—or on clearly defined spines attached to dendritic shafts. The remaining synapses were classified as targeting non-spine elements. Comparisons of synaptic densities across groups were performed using the Wilcoxon rank-sum test. The distribution of spine versus non-spine synapses was analysed using two-sided Fisher’s exact tests. Statistical significance was defined as P < 0.05.

Expansion microscopy analyses

We adapted a recent approach leveraging expansion microscopy and combinatorial antigen barcoding²⁶. Seven-month-old organoids cultured in either CDM4 or APM were transduced with a mix of nine AAVs (8 × 10⁸ viral genomes each), each expressing a CAG-driven spaghetti-monster fluorescent protein²⁷ tagged with distinct epitopes (V5, MYC, Flag, Ollas, E, S, HSV, protein C or Ty1). After 3 weeks, organoids were fixed in 4% PFA, embedded in gelatin, cryosectioned at 50 µm and mounted on HistoBond slides. The sections were expanded using the Magnify protocol²⁸, including methacrolein anchoring, overnight gelation at 37 °C, sample homogenization and expansion through water washes. Fully expanded samples were re-embedded in a non-expandable hydrogel to allow iterative immunostaining. Re-embedded samples were blocked and stained with primary antibodies overnight, washed and incubated with secondary antibodies the following night. For multiround immunostaining, antibody destaining was performed under the same homogenization conditions used for the initial sample preparation. Before imaging, samples were transferred to six-well glass-bottom plates. A Nikon W1-SORA spinning-disk confocal microscope was used: ×4 overview images identifying dense regions were followed by ×10 z stacks and ×40 tile scans. Large ×40 scans were stitched using BigStitcher; registration across rounds was performed with BigStream using a two-stage alignment pipeline (coarse RANSAC-based feature matching, followed by fine intensity-based registration using mutual maximum information). Affine transformation matrices were applied to full-resolution datasets. Preprocessing for low-signal-to-noise-ratio rounds included percentile-based contrast enhancement and 3D median filtering. Aligned volumes were converted to .wkw format, imported into Webknossos and merged into a single nine-channel dataset. Ten neurons per sample (CDM4, APM) were manually skeletonized based on multitag expression in the cytosol. Resulting traces (30 per condition) were converted to .swc files, and morphometric parameters were quantified using Navis in Python. In the Strahler analysis of CDM4 and APM neurons, terminal branches receive Strahler number one. When two branches with the same Strahler number meet at a branch point, the Strahler number of the next segment increases by one.

Dissociation of brain organoids and scRNA-seq

Organoids were dissociated as previously described⁹. Dissociated cells were resuspended in 100 μl of PBS, filtered with a 35 μm cell strainer tube (Corning) to remove aggregates and counted with a Countess II automated haematocytometer (Thermo Fisher Scientific) or a Cellometer K2 (Nexcelom Bioscience) and AOPI stain (Nexcelom Bioscience, CS2-0106-5ML). Single-cell suspensions were loaded into the Chromium Next GEM Chip G (10x Genomics, 1000120) and run with the Chromium Controller to generate single-cell GEMs. scRNA-seq libraries were prepared with the Chromium Single Cell 3′ Library and Gel Bead Kit v3.1 (10x Genomics, 1000268). The resulting libraries were pooled based on molar concentration and sequenced on a NovaSeq X or NovaSeq 6000 instrument (Illumina) either 28 bases for read 1, 55 bases for read 2 and 8 bases for index read 1, or 28 bases for read 1, 90 bases for read 2 and 10 bases for index read 1 (Supplementary Table 1). If necessary, after the first round of sequencing, we repooled libraries based on the actual number of cells in each library, and resequenced with the goal of producing an equal number of reads per cell for each sample (with a target read depth of 20,000 reads per cell).

scRNA-seq data analysis

The preprocessing of all scRNA-seq samples, including the conversion of raw BCL files into FASTQ files, reads alignment, and the generation of gene-barcode matrices, was performed using 10x Genomics Cell Ranger suite tool⁶⁴ following the steps described in our previous publication⁹. For de novo generated samples with moderate to high ambient RNA levels, CellBender⁶⁵ v.0.3.0 was used to remove systematic biases and background noises (Supplementary Table 1). The remove-background pipeline with the default parameters was used by taking the raw gene-by-cell count matrices produced by Cell Ranger pipelines as inputs. The filtered gene–barcode count matrices generated by Cell Ranger and CellBender were imported into Seurat (v.4.3.0)⁶⁶ for downstream analysis.

Genetic demultiplexing

In the cases in which multiple samples of different genetic backgrounds were pooled together and sequenced in the same lane to maximize the sequencing capacity, Demuxlet (v.1.0)⁶⁷ was applied with default parameters to determine the sample identity for each barcoded droplet and identify ambiguous droplets and doublets that contain two heterogeneous cells. Ambiguous droplets and doublets identified were removed from downstream analysis. To run Demuxlet, the BAM file output from Cell Ranger, the filtered barcode list generated by the initial preprocessing pipeline and reference variant-call-format (VCF) files for each expected cell line were used as inputs. The reference VCF files were prepared according to the method described in our previous publication⁹.

To address doublet overestimation when applying Demuxlet to multiplexed samples with ambient RNAs, after using CellBender, Souporcell (v.2.5)⁶⁸ was suggested by the developer and was used with default parameters as an alternative method, which allows genotype-free demultiplexing and clusters cells by the genetic variants detected within the input reads. The same set of input data plus the GRCh38 human reference genome FASTA file were used. We followed recommendations to include the reference VCF that has genetic variants of pooled individuals to enhance accuracy. Lastly, the Assign_Indiv_by_Geno.R function from the Demuxafy⁶⁹ framework was used to correlate clusters with donors in the reference SNP genotypes.

scRNA-seq data quality control, normalization, dimensionality reduction and clustering

For each scRNA-seq dataset, we performed initial quality-control filtering and retained cell profiles with total UMI counts of greater than 500 and less than 20,000, unique feature counts greater than 200 and a percentage of mitochondrial RNA less than 15. Seurat’s SCTransform function was used to normalize the data and regress out the mitochondrial RNA percentage. The RunPCA function was used to perform principal components analysis (PCA), and the top 30 principal components were used by the FindNeighbors function to construct a k-nearest-neighbour graph with the default parameters. Cells were then clustered by the Louvain algorithm using the FindClusters function with resolution = 0.3 for a rough granularity to start with. Clusters showing atypically low UMI counts, low unique feature counts and a high percentage of mitochondrial RNA were labelled as low-quality cells and removed (Supplementary Table 5).

Data compilation of group-published data

To unify data processing across datasets, raw scRNA-seq FASTQ files for organoids from refs. ^6,7 were reprocessed using Cell Ranger v.7.2.0 (the most notable change from the original processing is the inclusion of intronic reads in the count matrices).

Cluster annotation and data integration

Before any data merging/data integration, for each scRNA-seq data, we manually annotated each cluster by an ensemble approach, as described in our previous publication⁹. Doing so provides a solid way for cell type assignments by taking the previous domain knowledge, Seurat’s FindMarkers output and alignment with notations in previously published reference maps⁶ into consideration.

In pursuit of comparing organoids cultured under the conventional CDM4 and the novel APM conditions, we profiled APM-treated organoids at 4 months, 6 months, 9 months, 1 year and 1.5 years in this study. To facilitate downstream comparative analysis and enable more interpretable results, we refined cell type annotations between treatment conditions and corrected for sequencing batches as follows: For each age, scRNA-seq datasets from two conditions were merged by using Seurat’s merge function. Considering the notable data imbalance between CDM4- and APM-treated organoids at 4 months and to control the overall cell number disparity between conditions, a 50% downsampling was applied on the CDM4 scRNA-seq dataset while retaining the original sample proportions before merging. Features that are repeatedly variable across datasets in the list were obtained by using Seurat’s SelectIntegrationFeatures using the default parameters and were set as the variable features of the merged object. The RunPCA function was applied, and the merged data were then batch-corrected using Harmony⁷⁰ v.1.2.0. The downstream clustering used the standard Seurat pipeline by taking Harmony-corrected embeddings with a resolution set between 0.8 and 1.0. To examine the top upregulated genes within each cluster for cell type annotation, before using FindAllMarkers, Seurat’s PrepSCTFindMarkers was used to recalculate the SCT assay by taking the minimum median UMI across all datasets being merged as a scale factor to correct for sequencing depth variations across different datasets. To achieve higher resolution and finer granularity in certain clusters, including clusters that were labelled as cycling cells and clusters with underlying heterogeneity of cell types, subclustering was performed (Supplementary Table 5). Although the neuronal cluster identified at 5 years contained both cells expressing inhibitory markers and cells expressing excitatory markers, due to the small size of this cluster (101 cells), we did not attempt further subclustering. Owing to the absence of organoid samples with matched genetic backgrounds across conditions, the integrated 4-month data were excluded from the age-wise comparative analysis across conditions. Thereby, we retained cell profiles from organoids with matched genetic backgrounds for other timepoints for data representation and downstream analysis as follows: 4 months (CDM4: n = 6, 15,972 cells; APM: n = 2, 4,579 cells), 6 months (CDM4: n = 10, 64,817 cells; APM: n = 9, 41,231 cells), 9 months (CDM4: n = 6, 18,104 cells; APM: n = 3, 8,533 cells), 12 months (CDM4: n = 4, 7,802 cells; APM: n = 3, 6,084 cells), 1.5 years (CDM4: n = 2, 3,952 cells; APM: n = 2, 8,670 cells).

To assess the reproducibility of organoids, we calculated Aitchison distances based on cell type compositional data between every unique pair of samples at each timepoint by using the aDist function from the R package robCompositions (v.2.3.1)⁷¹. Aitchison distances between replicates within the treatment protocol were compared by two-sided Wilcoxon rank-sum tests. Five previously published human fetal cortex datasets^6,38,39, as noted in our previous publication⁹, were used to provide referential median sample-to-sample variabilities by their respective cell type annotations. Samples from the combined perinatal datasets from the Kriegstein lab^10,11 were grouped by age range (fetal samples by gestational month, and postnatal samples into 0–1 year, 1–2 years and 2–4 years) so that the Aitchison distances were calculated only between samples at similar ages.

To obtain a comprehensive overview of the organoid development across timepoints, we generated a developmental time map by merging the annotated scRNA-seq data of CDM4-treated organoids, including 424,720 cells from 110 organoid samples across 15 days (n = 1; 8,831 cells), 23 days (n = 6; 29,769 cells), 1 month (n = 6; 37,876 cells), 1.5 months (n = 6; 35,898 cells), 2 months (n = 6; 37,194 cells), 3 months (n = 31; 69,719 cells), 4 months (n = 6; 33,825 cells), 5 months (n = 6; 32,115 cells), 6 months (n = 12; 74,172 cells), 9 months (n = 10; 32,550 cells), 12 months (n = 4; 7,802 cells), 1.5 years (n = 6; 8,515 cells), 24 months (n = 4; 11,228 cells), 36 months (n = 2; 2,481 cells), 42 months (n = 2; 691 cells) and 60 months (n = 2; 2,054 cells). In general, we annotated 23 cell types, and we grouped cell types into several populations for broader representation as follows: (1) progenitors (aRG, oRGs, IPs); (2) excitatory neurons (ExNeu) (CFuPNs, CPNs, unspecified projection neurons (PNs) and newborn deep-layer projection neurons (newborn DL PNs)); (3) interneurons (INT) (immature interneurons, interneuron progenitors); (4) glial progenitors (GlialProg) (glial precursors, tRGs); (5) astrocytes (AST); and (6) OPCs. All other annotated cell types were grouped as ‘other’. Annotation of de novo generated datasets was performed manually, guided by the most current biological evidence from the literature and informed by our previously published findings. For the oldest organoid timepoints, where clusters were not always readily distinguishable based on canonical marker expression alone, annotations were assigned using an integrative approach that incorporated marker profiles, label transfer from earlier timepoints and biological judgment informed by the literature.

The merged data were jointly renormalized by SCTransform and reprocessed using the standard Seurat pipeline. The time-series object of APM-treated organoids from 4 months to 1.5 years was generated by the same method after merging the APM counterpart from age-wise integrated data as described above.

Analysis of previously published human fetal data

A combination of two previously published perinatal single-nucleus RNA-seq (snRNA-seq) datasets for human cortical development were used for comparison. The dataset from ref. ¹⁰ includes profiles from 709,372 nuclei of 169 brain tissue samples from 106 individuals, with a wide age range spanning from the second trimester of gestation to adulthood. The snRNA-seq data and metadata were downloaded from https://cells.ucsc.edu/?ds=pre-postnatal-cortex. We subsetted the snRNA-seq data to an age range from the second trimester up to 4 years, which is a proxy for the time window in our time-series data. The data were further subset by brain regions including cortex and frontal cortex with a spectrum of cell type populations including progenitors, excitatory neurons, interneurons, glial progenitors, astrocytes and OPCs. The dataset was further downsampled to include a maximum of 10,000 nuclei from any given annotated cell type (as defined by the original authors) at any given age range (second trimester, third trimester, 0–1 year, 1–2 years or 2–4 years), yielding a final subset of 151,749 nuclei from 56 individuals. Samples from ref. ¹¹ that had an estimated age of less than 100 post-conceptional days and came from cortex or prefrontal cortex regions were selected to be combined with the ref. ¹⁰ data. In total, 36,751 cells from 6 individuals met these criteria, and were merged with the above dataset. The combined reference, comprising 188,500 cells from 62 individuals, was collectively renormalized using Seurat’s SCTransform function. To compare with the reference human perinatal snRNA-seq dataset and explore what age the time-series organoids data maps to, we used labels by concatenating the information of cell type classes and age and performed reference-based label transfer using Seurat’s FindTransferAnchors and MapQuery functions.

Identification of MCPs associated with ageing

To explore the maturation profiles of the time-series organoids data, we used the aforementioned human perinatal dataset as inputs and applied DIALOGUE (v.1.0)¹², a tool designed to identify MCPs of co-regulated genes across different cell types and map how the transcriptome varies with changes in its environment. Following DIALOGUE’s vignette, a list of cell type objects was generated by the make.cell.type function to represent major cell populations in the human perinatal dataset by using a series of input arguments including the single-cell gene expression data, sample identity and PCA embeddings for each cell type of interest. Metadata of the SCTransform-normalized total UMI counts and the developmental age in days were included as a technical variability and a meaningful biological covariate to adjust for the analysis. Next, the DIALOGUE.run function was applied to the above list. As a result, five MCPs were identified, and each MCP output cell-type-specific upregulated and downregulated gene sets. We ran Seurat’s AddModuleScores function on the combined perinatal tissue reference, using each of DIALOGUE’s MCP and cell-type-specific gene lists as inputs. For each MCP, we assigned an overall score to each cell by subtracting the module score calculated from the downregulated gene set that matched that cell’s cell type from the score of the corresponding upregulated gene set. In perinatal tissues, of the five MCPs, MCP4 most consistently showed the highest correlation with age within most cell population and across all cells; we therefore focused on the cell-type-specific up-minus-down MCP4 module scores as maturation scores and we adopted the same approach to calculate the maturation scores by using the above master gene sets in the context of our 15-day to 5-year organoid timeseries dataset, our 4-month to 18-month APM versus CDM4 organoid dataset, and our human heterochronic and monochronic organoid datasets.

For scatter plot representations, the R geom_smooth() function was used to add smoothed conditional means (black lines) with default 95% CIs (grey regions). Points were given slight x-axis jitter to improve visibility; horizontal position within a timepoint does not indicate differences in age. Heat map representations were generated to show the expression of either all MCP4 genes (in organoids) or cell-type-specific MCP4 genes (within a given cell population in organoids), respectively, over time. Values for each timepoint represent the mean of the counts-per-million of the pseudobulked count matrices for either each organoid or a given cell type within each organoid of that timepoint, respectively. z-Score normalization was performed for each gene. Genes are ordered by the weighted expression over time, by positioning genes expressed mostly at earlier time points at the top, and genes expressed mostly at later time points at the bottom.

GSEA

To investigate and compare temporal changes of astrocytes in organoids and temporal changes of astrocytes in endogenous perinatal human cortex, differential expression analyses were first performed in each dataset, treating age as a continuous variable. All genes were ranked based on the significance and direction of change (using the negative-log-transformed P value multiplied by the sign of the log-transformed fold change as the ranking metric). The top 500 upregulated genes over time in perinatal tissue astrocytes and the top 500 downregulated genes over time (marked with vertical blue lines) were used as gene sets for Kolmogorov–Smirnov tests based on the ranking in organoids. For visualizing results, the running enrichment scores (y axis) are shown as red (for upregulated genes) and blue (for downregulated genes) curves, which track the cumulative walk down the ranked data by increasing/decreasing given the presence/absence of a gene in the target gene set. Maximum absolute values of the running enrichment scores are the Kolmogorov–Smirnov tests’ D statistic. P values of the Kolmogorov–Smirnov tests are shown in boxed text. Identical analyses were performed by taking the endogenous tissue data as the rank data, and the top 500 up and downregulated genes as gene sets.

Analysis of presynapse, synapse and post-synapse signature modules

Gene sets of GO terms: synapse (GO: 0045202), presynapse (GO: 0098793) and post-synapse (GO: 0098794) containing 1,778, 879 and 1,181 genes, respectively, were retrieved from the SynGO knowledge base of release v.1.2 (https://www.syngoportal.org/)⁷². Each gene set was used to calculate module scores for cells in integrated scRNA-seq data of CDM4- and APM-treated organoids at 6 months, 9 months and 12 months using Seurat’s AddModuleScore function. To compare, under each age, whether organoids under two treatment conditions displayed different averaged expression levels for each program, LMEMs were built by the R package lme4⁷³ to model the scores by including which organoid the cell belongs to and the genetic background of that organoid as random effects (score ~ (1|organoid) + (1|genotype)). To adjust for cell count differences in each broad cell population across organoids, the reciprocal of the square root of the observed cell count per population per organoid was included as a weight factor. Models with or without treatment as an additive fixed effect were compared by the anova function from the R stats package and P values across tests of the above three signature modules were adjusted using the Benjamini–Hochberg method.

Analysis of Hallmark apoptosis, hypoxia and glycolysis gene signatures

We used the R package msigdbr v.25.1.0 to access MSigDB Hallmark gene sets for apoptosis, hypoxia and glycolysis. Each of these gene sets was used as an input to Seurat’s AddModuleScores function to assign scores to cells from our 15-day to 5-year organoid timeseries dataset, as well as to our 4-month to 18-month APM versus CDM4 organoid dataset. For each Hallmark module score, we used mixed-effects models through lme4 to assess whether there were significant changes over time, both globally (across all cell types) by comparing the model score ~ log(months) + (1|organoid) + (1|broad_cell_type) against a similar model but excluding a covariate for the age of the organoids. Similar tests were performed within each broad cell type independently (excluding the covariate for cell type, of course).

Comparing cell type proportions between treatment conditions

The cell type proportion changes between CDM4- and APM-treated organoids at 4 months, 6 months, 9 months, 1 year and 1.5 years were evaluated by NBME models. A matrix of cell counts for each annotated cell type per organoid was generated. Given a cell type of interest, a NBME model was created using the glmer.nb function from lme4⁷³, which takes the genetic background of an organoid replicate as a random effect to account for the variability in genotypes. The model also includes the logarithm of the total cell number of each organoid sample as an offset to account for the cell count differences (model formula: cell type ~ offset(log(library size)) + (1|genotype)). The significance of cell type proportion differences between organoids under two treatment conditions was calculated by the one-sided likelihood ratio test comparing negative binomial models with and without treatment as an additive fixed effect using the anova function. The Benjamini–Hochberg method was used to adjust P values from multiple independent tests. Cell type proportion differences between MonChr chimeroids and 9 month organoids were calculated using the same methods as above and only considering cell types with a minimum abundance of 5% in 9 month organoids.

Comparing cell cycle proportions between treatment conditions

Following Seurat’s vignette, we used the CellCycleScoring function to calculate G2M and S-phase module scores per cell in our APM versus CDM4 organoid dataset. Cells that scored below 0.1 in both G2M and S-phase modules were designate as non-cycling, and cells that scored 0.1 or higher in either module were designated as cycling. To test differences in the proportion of progenitors that are cycling between APM-treated and CDM4-treated organoids, we considered each type of progenitor cell independently, and we constructed binomial mixed-effects models vi lme4 with the design: cycling ~ months + treatment + (1|organoid), compared against a similar model excluding the covariate for treatment using the anova function.

DEG analysis and GO analysis

To identify differentially expressed genes (DEGs) between CDM4- and APM-treated organoids at each time point (6, 9 and 12 months), we used a pseudo-bulk approach to analyse the data at the sample level. We took pan-ExNeu, AST and IN as partitions of interest by using the above-mentioned grouping criteria. Pseudobulk profiles were generated on a partition basis by aggregating UMI counts for each gene across cells of that partition within each organoid sample. We require at least 20 cells within a partitioned cluster per organoid to retain the sample. Genes with fewer than ten total UMI counts in at least two samples were excluded. The pseudo-bulk differential expression analysis was performed by using DESeq2⁷⁴ with the design formula ~genotype + treatment to identify DEGs between treatment conditions while regressing out the variation by genetic backgrounds. DEGs were defined by a cut-off with the absolute value of log₂-transformed fold change > 0.5 and FDR-adjusted P < 0.05. Upregulated DEGs were used by the enrichGO function from clusterProfiler⁷⁵ to perform GO enrichment analysis. The simplify function with a default similarity cutoff of 0.7 from the same package was used to remove the semantic redundancy of enriched GO terms.

Comparative analysis of DIALOGUE-determined maturation scores

To examine whether APM-treated organoids show different DIALOGUE maturation scores in a time-series or time-specific manner, the age-wise integrated data from 4 months to 1 year were merged and jointly renormalized by SCTransform and reprocessed by the standard Seurat pipeline. We used Seurat’s AddModuleScore function to calculate module scores using the DIALOGUE MCP4 full list of ‘up’ and ‘down’ genes for each broad cell type, separately, and then took the subtraction of the corresponding cell type’s ‘up − down’ for each cell as a single maturation score, as described above. For the comparative analysis, the z score, with a mean of 0, was used to represent the magnitude. Maturation module scores across the entire dataset were standardized using the R function scale. Focusing on the excitatory neurons, we first assessed the treatment effect over time by treating age as a continuous variable. The z scores were averaged across all cells within the partition of interest per organoid. A baseline linear model was created to evaluate how maturation scores change with age (score ~ age). A second linear model incorporating treatment in addition (score ~ treatment + age) was used to compare with the baseline model by anova function from the R stats package to test whether adding treatment as an additive effect significantly improves the model fitting. Next, to assess in a time-specific manner, LMEMs were used to model the scores by including which organoid the cell belongs to and the genetic background of that organoid as random effects (score ~ (1|organoid) + (1|genotype)) and a weight factor to adjust for cell count differences in excitatory neurons across organoids. Models with or without treatment as an additive fixed effect were evaluated by anova function. The above statistical models were built by the R lme4⁷³ package.

Confidence intervals

Throughout the text, 95% CIs for statistical tests were calculated using the confint function from the R stats package v.4.4.0, except in the case of the results from Fisher’s exact tests, for which confidence intervals are reported directly from the outputs of the fisher.test function from the stats package.

Reporting summary

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

Code used during data analysis is available at GitHub (https://github.com/tfaits/Arlotta_Lab_LongTermOrganoids), with the exception of custom code used for MEA analysis, which is available on request.

1. Ciceri, G. et al. An epigenetic barrier sets the timing of human neuronal maturation. Nature 626, 881–890 (2024).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

2. Del Dosso, A., Urenda, J.-P., Nguyen, T. & Quadrato, G. Upgrading the physiological relevance of human brain organoids. Neuron 107, 1014–1028 (2020).

   Article  PubMed  PubMed Central  Google Scholar 

3. Casimir, P., Iwata, R. & Vanderhaeghen, P. Linking mitochondria metabolism, developmental timing, and human brain evolution. Curr. Opin. Genet. Dev. 86, 102182 (2024).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

4. Gordon, A. et al. Long term maturation of human cortical organoids matches key early postnatal transitions. Nat. Neurosci. 24, 331–342 (2021).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

5. Fenlon, L. R. Timing as a mechanism of development and evolution in the cerebral cortex. Brain Behav. Evol. 97, 8–32 (2021).

   Article  PubMed  Google Scholar 

6. Uzquiano, A. et al. Proper acquisition of cell class identity in organoids allows definition of fate specification programs of the human cerebral cortex. Cell 185, 3770–3788 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

7. Velasco, S. et al. Individual brain organoids reproducibly form cell diversity of the human cerebral cortex. Nature 570, 523–527 (2019).

   Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

8. He, Z. et al. An integrated transcriptomic cell atlas of human neural organoids. Nature 635, 690–698 (2024).

9. Antón-Bolaños, N. et al. Brain Chimeroids reveal individual susceptibility to neurotoxic triggers. Nature 631, 142–149 (2024).

   Article  ADS  PubMed  PubMed Central  Google Scholar 

10. 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 

11. Wang, L. et al. Molecular and cellular dynamics of the developing human neocortex. Nature 647, 169–178 (2025).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

12. Jerby-Arnon, L. & Regev, A. DIALOGUE maps multicellular programs in tissue from single-cell or spatial transcriptomics data. Nat. Biotechnol. 40, 1467–1477 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

13. Moqri, M., Poganik, J. R., Horvath, S. & Gladyshev, V. N. What makes biological age epigenetic clocks tick. Nat. Aging 5, 335–336 (2025).

    Article  PubMed  PubMed Central  Google Scholar 

14. Price, A. J. et al. Divergent neuronal DNA methylation patterns across human cortical development reveal critical periods and a unique role of CpH methylation. Genome Biol. 20, 196 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

15. Luo, C. et al. Cerebral organoids recapitulate epigenomic signatures of the human fetal brain. Cell Rep. 17, 3369–3384 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

16. Christensen, B. C. et al. Aging and environmental exposures alter tissue-specific DNA methylation dependent upon CpG island context. PLoS Genet. 5, e1000602 (2009).

    Article  PubMed  PubMed Central  Google Scholar 

17. Horvath, S. DNA methylation age of human tissues and cell types. Genome Biol. 14, R115 (2013).

    Article  PubMed  PubMed Central  Google Scholar 

18. Steg, L. C. et al. Novel epigenetic clock for fetal brain development predicts prenatal age for cellular stem cell models and derived neurons. Mol. Brain 14, 98 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

19. Bardy, C. et al. Neuronal medium that supports basic synaptic functions and activity of human neurons in vitro. Proc. Natl Acad. Sci. USA 112, E2725–2734 (2015).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

20. Perriot, S., Canales, M., Mathias, A. & Du Pasquier, R. Differentiation of functional astrocytes from human-induced pluripotent stem cells in chemically defined media. STAR Protoc. 2, 100902 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

21. Pavarino, E. C. et al. mEMbrain: an interactive deep learning MATLAB tool for connectomic segmentation on commodity desktops. Front. Neural Circuits 17, 952921 (2023).

22. Dhainaut, M. et al. Spatial CRISPR genomics identifies regulators of the tumor microenvironment. Cell 185, 1223–1239 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Kudo, T., Lane, K. & Covert, M. W. A multiplexed epitope barcoding strategy that enables dynamic cellular phenotypic screens. Cell Syst. 13, 376–387 (2022).

    Article  CAS  PubMed  Google Scholar 

24. Rovira-Clavé, X. et al. Spatial epitope barcoding reveals clonal tumor patch behaviors. Cancer Cell 40, 1423–1439 (2022).

    Article  PubMed  PubMed Central  Google Scholar 

25. Wroblewska, A. et al. Protein barcodes enable high-dimensional single-cell CRISPR screens. Cell 175, 1141–1155 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

26. Park, S. Y. et al. Combinatorial protein barcodes enable self-correcting neuron tracing with nanoscale molecular context. Preprint at bioRxiv https://doi.org/10.1101/2025.09.26.678648 (2025).

27. Viswanathan, S. et al. High-performance probes for light and electron microscopy. Nat. Methods 12, 568–576 (2015).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

28. Klimas, A. et al. Magnify is a universal molecular anchoring strategy for expansion microscopy. Nat. Biotechnol. 41, 858–869 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

29. Berg, J. et al. Human neocortical expansion involves glutamatergic neuron diversification. Nature 598, 151–158 (2021).

30. Paulsen, B. et al. Autism genes converge on asynchronous development of shared neuron classes. Nature 602, 268–273 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

31. Pachitariu, M., Sridhar, S., Pennington, J. & Stringer, C. Spike sorting with Kilosort4. Nat. Methods 21, 914–921 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

32. Oberst, P., Agirman, G. & Jabaudon, D. Principles of progenitor temporal patterning in the developing invertebrate and vertebrate nervous system. Curr. Opin. Neurobiol. 56, 185–193 (2019).

    Article  CAS  PubMed  Google Scholar 

33. Shen, Q. et al. The timing of cortical neurogenesis is encoded within lineages of individual progenitor cells. Nat. Neurosci. 9, 743–751 (2006).

    Article  CAS  PubMed  Google Scholar 

34. Mathiesen, S. N., Lock, J. L., Schoderboeck, L., Abraham, W. C. & Hughes, S. M. CNS transduction benefits of AAV-PHP.eB over AAV9 are dependent on administration route and mouse strain. Mol. Ther. Methods Clin. Dev. 19, 447–458 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

35. Fortuna, M. G. et al. AAV-PHP.eB achieves superior neuronal transduction over AAV9 in pigtail macaques following intracerebroventricular administration. Mol. Ther. Methods Clin. Dev. 33, 101636 (2025).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

36. Chan, K. Y. et al. Engineered AAVs for efficient noninvasive gene delivery to the central and peripheral nervous systems. Nat. Neurosci. 20, 1172–1179 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

37. Mitchell, J. M. et al. Mapping genetic effects on cellular phenotypes with “cell villages”. Preprint at bioRxiv https://doi.org/10.1101/2020.06.29.174383 (2020).

38. Trevino, A. E. et al. Chromatin and gene-regulatory dynamics of the developing human cerebral cortex at single-cell resolution. Cell 184, 5053–5069 (2021).

    Article  CAS  PubMed  Google Scholar 

39. Polioudakis, D. et al. A single-cell transcriptomic atlas of human neocortical development during mid-gestation. Neuron 103, 785–801 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

40. Quadrato, G. et al. Cell diversity and network dynamics in photosensitive human brain organoids. Nature 545, 48–53 (2017).

41. Grosswendt, S. et al. Epigenetic regulator function through mouse gastrulation. Nature 584, 102–108 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

42. Sampath Kumar, A. et al. Spatiotemporal transcriptomic maps of whole mouse embryos at the onset of organogenesis. Nat. Genet. 55, 1176–1185 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

43. Frey, E. C., Humm, J. L. & Ljungberg, M. Accuracy and precision of radioactivity quantification in nuclear medicine images. Semin. Nucl. Med. 42, 208–218 (2012).

    Article  PubMed  PubMed Central  Google Scholar 

44. Haggerty, C. et al. Dnmt1 has de novo activity targeted to transposable elements. Nat. Struct. Mol. Biol. 28, 594–603 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

45. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12 (2011).

    Article  Google Scholar 

46. Xi, Y. & Li, W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinform. 10, 232 (2009).

    Article  Google Scholar 

47. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).

    Article  PubMed  PubMed Central  Google Scholar 

48. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

49. Sun, D. et al. MOABS: model based analysis of bisulfite sequencing data. Genome Biol. 15, R38 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

50. Jühling, F. et al. metilene: fast and sensitive calling of differentially methylated regions from bisulfite sequencing data. Genome Res. 26, 256–262 (2016).

    Article  PubMed  Google Scholar 

51. McLean, C. Y. et al. GREAT improves functional interpretation of cis-regulatory regions. Nat. Biotechnol. 28, 495–501 (2010).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

52. Charlton, J. et al. TETs compete with DNMT3 activity in pluripotent cells at thousands of methylated somatic enhancers. Nat. Genet. 52, 819–827 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

53. Hetzel, S. et al. Acute lymphoblastic leukemia displays a distinct highly methylated genome. Nat. Cancer 3, 768–782 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

54. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

55. Zhou, W. et al. DNA methylation loss in late-replicating domains is linked to mitotic cell division. Nat. Genet. 50, 591–602 (2018).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

56. Pelegí-Sisó, D., de Prado, P., Ronkainen, J., Bustamante, M. & González, J. R. methylclock: a Bioconductor package to estimate DNA methylation age. Bioinformatics 37, 1759–1760 (2021).

    Article  PubMed  Google Scholar 

57. Shireby, G. L. et al. Recalibrating the epigenetic clock: implications for assessing biological age in the human cortex. Brain 143, 3763–3775 (2020).

    Article  PubMed  PubMed Central  Google Scholar 

58. Gu, Z., Eils, R., Schlesner, M. & Ishaque, N. EnrichedHeatmap: an R/Bioconductor package for comprehensive visualization of genomic signal associations. BMC Genom. 19, 234 (2018).

    Article  Google Scholar 

59. Kirillov, A. NeuroExplorer Manual (Nex Technologies, 1998–2014).

60. Kasthuri, N. et al. Saturated reconstruction of a volume of neocortex. Cell 162, 648–661 (2015).

61. Hayworth, K. J. et al. Imaging ATUM ultrathin section libraries with WaferMapper: a multi-scale approach to EM reconstruction of neural circuits. Front. Neural Circuits 8, 68 (2014).

62. Berger, D. R., Seung, H. S. & Lichtman, J. W. VAST (Volume Annotation and Segmentation Tool): efficient manual and semi-automatic labeling of large 3D image stacks. Front. Neural Circuits 12, 88 (2018).

63. Ronneberger, O., Fischer, P. & Brox, T. U-Net: Convolutional networks for biomedical image segmentation. In Proc. Medical Image Computing and Computer-Assisted Intervention 2015 (eds Navab, N. et al.) 234–241 https://doi.org/10.1007/978-3-319-24574-4_28 (Springer, 2015).

64. Zheng, G. X. Y. et al. Massively parallel digital transcriptional profiling of single cells. Nat. Commun. 8, 14049 (2017).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

65. 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 

66. Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

67. Kang, H. M. et al. Multiplexed droplet single-cell RNA-sequencing using natural genetic variation. Nat. Biotechnol. 36, 89–94 (2018).

    Article  CAS  PubMed  Google Scholar 

68. Heaton, H. et al. Souporcell: robust clustering of single-cell RNA-seq data by genotype without reference genotypes. Nat. Methods 17, 615–620 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

69. Neavin, D. et al. Demuxafy: improvement in droplet assignment by integrating multiple single-cell demultiplexing and doublet detection methods. Genome Biol. 25, 94 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

70. 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 

71. Templ, M., Hron, K. & Filzmoser, P. in Compositional Data Analysis 341–355 https://doi.org/10.1002/9781119976462.ch25 (John Wiley & Sons, 2011).

72. Koopmans, F. et al. SynGO: an evidence-based, expert-curated knowledge base for the synapse. Neuron 103, 217–234 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

73. Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).

    Article  Google Scholar 

74. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

75. Yu, G., Wang, L.-G., Han, Y. & He, Q.-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Author notes

1. Irene Faravelli

   Present address: Department for Clinical Sciences and Community Health (DISCCO), University of Milan, Milan, Italy

2. Noelia Antón-Bolaños

   Present address: Department of Translational Neuroscience, University Medical Center Utrecht (UMCU), Brain Center, Utrecht University, Utrecht, The Netherlands

3. These authors contributed equally: Irene Faravelli, Noelia Antón-Bolaños

Authors and Affiliations

1. Department of Stem Cell & Regenerative Biology, Harvard University, Cambridge, MA, USA

   Irene Faravelli, Noelia Antón-Bolaños, Anqi Wei, Tyler Faits, Abhishek Sampath Kumar, Sophia Andreadis, Rahel Kastli, Marta Montero Crespo, Sayara Silwal, Sung Min Yang & Paola Arlotta

2. Stanley Center for Psychiatric Research, Broad Institute of MIT and Harvard, Cambridge, MA, USA

   Irene Faravelli, Noelia Antón-Bolaños, Anqi Wei, Tyler Faits, Abhishek Sampath Kumar, Sophia Andreadis, Rahel Kastli, Marta Montero Crespo, Sung Min Yang, Xian Adiconis, Joshua Z. Levin & Paola Arlotta

3. Department of Genome Regulation, Max Planck Institute for Molecular Genetics, Berlin, Germany

   Mara Steiger, Alexander Kovacsovics, Helene Kretzmer & Alexander Meissner

4. Digital Health Cluster, Hasso Plattner Institute for Digital Engineering, Digital Engineering Faculty, University of Potdsdam, Potsdam, Germany

   Mara Steiger & Helene Kretzmer

5. McGovern Institute, Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA

   Daniel Leible, Elizabeth Zhang, Bobae An & Edward S. Boyden

6. Department of Molecular and Cellular Biology and The Center for Brain Science, Harvard University, Cambridge, MA, USA

   Yaron Meirovitch & Jeff Lichtman

7. Klarman Cell Observatory, Broad Institute of MIT and Harvard, Cambridge, MA, USA

   Xian Adiconis, Joshua Z. Levin & Aviv Regev

8. Genentech, San Francisco, CA, USA

   Aviv Regev

9. Broad Institute of MIT and Harvard, Boston, MA, USA

   Alexander Meissner

Contributions

I.F. and N.A.-B. are equally contributing first authors; A.W. and T.F. are equally contributed second authors. N.A.-B., I.F. and P.A. conceived and designed the experiments. N.A.-B., I.F. and S.A. generated, cultured and characterized all of the organoids in this study, with the help of S.S.; P.A. supervised their work and contributed to data interpretation. N.A.-B. and I.F. worked on cell type assignments and scRNA-seq data analysis, and A.W. and T.F. performed all of the computational work under the supervision of A.R. and P.A.; A.S.K. performed methylome experiments and helped with analysis and interpretation. A.K. prepared the methylome libraries for sequencing. M.S. performed the computational analyses for the methylome data under the supervision of H.K. and A.M.; N.A.-B., I.F., S.A. and S.S. performed IHC and image acquisition. S.S. performed IHC and image acquisition and quantification related to SATB2 and FOS staining. A.S.K. performed NeuN IHC on long-term organoids. R.K. performed MEA recordings and S.M.Y. performed data analysis with the help of R.K.; M.M.C. performed all of the EM experiments and ground-truth analysis with the help of S.A. and Y.M.; Y.M. worked on the quantification algorithm for synapse quantification under the supervision of J.L.; N.A.-B. and I.F. performed scRNA-seq experiments with help from X.A. and supervision of P.A. and J.Z.L.; D.L. designed and performed the expansion microscopy epitope tag barcoding experiments and analyses, with help from E.Z., and B.A. contributed to early development of the epitope tag barcoding constructs, all under the supervision of E.S.B.; N.A.-B. and A.S.K. developed mouse endogenous chimeroids and performed the related experiments. N.A.-B. and I.F. performed scRNA-seq library preparation for the vast majority of the preparations included in this study. N.A.-B., I.F. and P.A. wrote the manuscript with contributions from all of the authors. All of the authors read and approved the final manuscript.

Corresponding author

Correspondence to Paola Arlotta.

Competing interests

P.A. is a scientific advisory board member at Foresite Labs, and is a co-founder of Vesalius and a co-founder and equity holder at Foresite Labs and at Avatar Bio. A.R. is a founder and equity holder of Celsius Therapeutics, an equity holder in Immunitas Therapeutics and, until 31 August 2020, was a scientific advisory board member of Syros Pharmaceuticals, Neogene Therapeutics, Asimov and Thermo Fisher Scientific. From 1 August 2020, A.R. has been an employee of Genentech and has equity in Roche. A.M. is a co-founder and scientific advisor of Harbinger Health, and is on the scientific advisory board of Zymo Research.

Peer review information

Nature thanks Colette Dehay, Andrew Jaffe, Selene Lomoio and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Reporting Summary (download PDF )

Peer Review File (download PDF )

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