Directed differentiation of human hindbrain neuroepithelial stem cells recapitulates cerebellar granule neurogenesis

Dysregulation of cerebellar granule neurogenesis underlies many human disease states, including movement disorders and medulloblastoma (Kawamura et al., 2021; Schüller et al., 2008; Vanner et al., 2014; Yang et al., 2008). Our understanding of mammalian cerebellar granule neuron (CGN) development largely depends on mouse studies, and it has long been assumed that this process is roughly equivalent in humans (Burgoyne and Cambray-Deakin, 1988). However, it has been shown that human and mouse granule lineage development differ temporally and anatomically (Aldinger et al., 2021; Haldipur et al., 2019). Specifically, the human rhombic lip (RL) persists for longer in development than its mouse counterpart and exists as ventricular zone (RL-VZ) and sub-ventricular zone (RL-SVZ) compartments. Therefore, cells underlying human disease states may not have direct equivalents in the mouse brain, preventing them from being faithfully modeled in a mouse system.

CGN differentiation begins when bone morphogenetic protein (BMP) molecules, which are secreted by fourth ventricle roof plate cells, act on adjacent rhombomere 1 (r1) neuroepithelium to induce RL formation (Alder et al., 1999). The RL is the germinal zone of cerebellar glutamatergic neurons and is characterized by ATOH1 and OLIG3 expression (Akazawa et al., 1995; Lowenstein et al., 2021). ATOH1-expressing RL progenitors tangentially migrate over the nascent cerebellum to form a superficial external granule layer (EGL) (Ben-Arie et al., 1997; Miale and Sidman, 1961). Cells populating the EGL are called granule neuron progenitors (GNPs). Purkinje neuron-secreted sonic hedgehog (SHH) is a GNP mitogen and expands the GNP population in the EGL (Wechsler-Reya and Scott, 1999). SHH signaling protects ATOH1 from proteasomal degradation (Forget et al., 2014), allowing maintenance of the primary cilium and transduction of SHH signaling (Chang et al., 2019). Notch signaling is also active within GNPs in the EGL and prevents precocious differentiation by maintaining a proliferative state (Solecki et al., 2001). As development proceeds, promoters of GLI1 and GLI2, the primary downstream effectors of SHH signaling are deacetylated, repressing the pathway and inducing CGN differentiation (Tiberi et al., 2014). Contemporaneously, differentiating GNPs begin to produce brain-derived neurotrophic factor (BDNF), which acts in an autocrine and paracrine manner to initiate and guide radial migration to the inner granule layer (IGL) of the cerebellum (Borghesani et al., 2002). Granule neurons in the IGL extend axons up through the Purkinje layer to the molecular layer of the cerebellum. Here, their axons bifurcate into parallel fibers that form glutamatergic synapses with Purkinje neurons, which are ultimately responsible for relaying cerebellar output to the motor cortex via the deep cerebellar nuclei (Allen and Tsukahara, 1974).

Differentiation of mouse and human embryonic stem (ES) cells to cerebellar granule neurons has been reported (Behesti et al., 2021; Erceg et al., 2010; Salero and Hatten, 2007; Su et al., 2006). However, these studies focused on characterizing end products of differentiation, rather than developmental intermediaries, from which disease states are more likely to arise (Kawamura et al., 2021; Selvadurai et al., 2020; Vanner et al., 2014). Furthermore, they relied on co-culture of differentiating cells with mouse cerebellar cells or with growth medium conditioned by mouse cells. In such cases, it is unknown which microenvironmental factors are influencing differentiation. Functional characterization of differentiated neurons is also lacking in these works. Finally, it remains unclear whether a mature neuronal state is achieved in some of these studies, owing to reported expression of immature markers, such as GLI1, in the final differentiation timepoint (Erceg et al., 2010).

Human hindbrain neuroepithelial stem (hbNES) cells have been isolated from the embryonic hindbrain (Tailor et al., 2013). hbNES cells are long-term self-renewing stem cells that maintain their regional and temporal identities in vitro and have been shown to be responsive to various signaling cues. Previous studies on pluripotent stem cell-derived NES cells (lt-NES) have demonstrated the developmental plasticity of these cells by differentiating them to various neuronal types in vitro, including motor neurons and dopaminergic neurons (Falk et al., 2012; Koch et al., 2009). However, the ability of hbNES cells to recapitulate hindbrain developmental processes remains unknown. Given these past results, we hypothesized that hbNES cells should possess the potential for differentiation towards the CGN lineage. If this could be achieved, a new genetically tractable model system could be created with applications in understanding human cerebellar development and pathology.

Here, we report efficient differentiation of CGNs from four independently derived hbNES cell lines. Using bulk RNA sequencing (RNA-seq), single cell Multiome sequencing (scMultiome-seq) and electrophysiology, we show that differentiating hbNES cells transit through RL and GNP stages before functionally maturing. We observed recapitulation of important human-specific features, such as the presence of RL-VZ and RL-SVZ identities, and the persistence of RL-like cells through to the end of differentiation. At our final day 56 differentiation timepoint, we obtained CGNs at a frequency of 86% and observed unipolar brush cells (UBCs). Although UBCs have been previously described in human pluripotent stem cell-derived cerebellar organoids (Silva et al., 2021), we show using transcriptional lineage reconstruction that these cells likely arise from a bipotent RL precursor, as is the case for mice (Englund et al., 2006). We further use our system to investigate the role of SHH during granule neuron differentiation. In addition to functioning as a mitogen, high SHH levels aberrantly promote GABAergic fate. Interestingly, we uncover a potential pro-differentiation role for SHH in the CGN lineage as spontaneous neuronal activity is diminished in its absence. Our work establishes and characterizes a human system of CGN differentiation that recapitulates key features of this process, and is well suited to modeling development and diseases such as medulloblastomas and movement disorders.

We cultured hbNES cell lines isolated from four different human embryos aged gestational weeks 5 to 6 of both sexes (Tailor et al., 2013) (Table S1A). Under adherent conditions, we observed all hbNES cell lines self-organized into neural rosettes characterized by radial symmetry. hbNES cells expressed the neuroepithelial transcription factor (TF) PLZF and tight-junction protein ZO1 on their apical surfaces (Fig. 1A). In their undifferentiated state, we confirmed that hbNES cells were not committed to the neural or glial lineage, indicated by absence of βIII-tubulin and GFAP staining, respectively (Fig. 1B). Furthermore, hbNES cells maintain expression of a repertoire of neuroepithelial genes and transcription factors, such as SOX1 and SOX2 (Fig. S1A,B). Together, these results suggest maintenance of neuroepithelial temporal identity in vitro. Importantly, hbNES cells maintained a regional identity restricted to r1, based on high expression of EN1 and EN2 (Davis and Joyner, 1988), and GBX2 (Joyner et al., 2000), along with low or absent expression of forebrain (FOXG1, OTX1 and OTX2) (Joyner et al., 2000; Martynoga et al., 2005) and posterior hindbrain (HOX genes) (Gavalas et al., 2003) markers (Fig. 1C), as measured by RNA-seq. To interrogate the heterogeneity of our hbNES cell population, we profiled the transcriptional and epigenetic landscapes of 8370 undifferentiated hbNES cells in a representative line (SAI5) using scMultiome-seq. Through nearest neighbor clustering, we identified three clusters within this population (Fig. 1D). We observed ubiquitous expression of classical neuroepithelial genes LIN28A and LIN28B (Herrlinger et al., 2019), PRTG (Wong et al., 2010), and SOX2 throughout these clusters (Fig. 1E-H), validating our previously made observations from bulk RNA-seq and immunocytochemistry (ICC). Among this undifferentiated population, we discovered a cluster of cells in which PAX3 and TFAP2B were significantly differentially expressed (P=3.25e⁻¹¹⁵ and 1.05e⁻¹²⁹, respectively) (Fig. 1I,J and Table S2C). We annotated this cluster as ‘Primed-NES’, based on the role of these two TFs in the differentiation of GABAergic cerebellar neurons (Zainolabidin et al., 2017). This observation corroborates a previous report that hbNES cells generate GABAergic interneurons when differentiated spontaneously (Tailor et al., 2013). To further interrogate this Primed-NES population at the epigenetic level, we performed differential TF motif activity analysis of the Primed-NES cluster relative to NES-1 and NES-2 clusters (Table S2K). In agreement with our previous observation that Primed-NES cells are inclined towards the GABAergic fate, we found that binding motifs for the GABAergic cerebellar lineage TFs TFAP2B (motif ID: MA0812.1, P=0) and PTF1A (motif ID: MA1619.1, P=5.13e⁻¹⁴⁸) were significantly more active in Primed-NES cells (Fig. 1K,L). However, to our surprise, Primed-NES cells also had significantly active cerebellar glutamatergic lineage TF motifs for ATOH1 (motif ID: MA1467.1, P=5.68e⁻¹⁵⁴) and ZIC1/ZIC2 (Aruga et al., 2002) (motif ID: 1628.1, P=9.83e⁻¹⁴³) (Fig. 1M,N). Therefore, we conclude that although Primed-NES cells may be transcriptionally primed for GABAergic differentiation, epigenetically they remain plastic, and potentially can be differentiated towards the glutamatergic CGN lineage. Taken together, these data demonstrate that hbNES cells retain their regional and temporal identities in vitro, and a subpopulation of them is primed for specific neuronal lineage differentiation.

Having verified the integrity of our starting hbNES cell population, we next sought to differentiate hbNES cells to CGNs to generate a model system with which to study developmentally relevant intermediate states. We modified previously reported differentiation protocols using embryonic stem (ES) cells and induced pluripotent stem (iPS) cells as their starting populations to start with hbNES cells instead (Fig. 2A) (Erceg et al., 2012, 2010; Salero and Hatten, 2007). First, we withdrew epidermal growth factor (EGF) and fibroblast growth factor (FGF) from the culture medium and added BMP6, BMP7 and GDF7 for 7 days to induce RL identity. At this timepoint, we added SHH and jagged 1 (JAG1), a Notch2 ligand, to promote GNP formation and proliferation. We also supplemented the growth medium with BDNF and neurotrophin 3 (NT3) from 7 days onward to induce neurite growth and maintain cell survival (Ghosh et al., 1994; Segal et al., 1995).

We first conducted quantitative reverse-transcriptase PCR (qRT-PCR) gene expression analysis of differentiating hbNES cells at the timepoints described in Fig. 2A. We measured the expression of 24 genes spanning the entire CGN lineage from hbNES cell to mature neuron. We conducted hierarchical clustering of relative expression values and noted timepoints tended to cluster together across the four different cell lines, indicating a robust response of hbNES cells to our differentiation protocol (Fig. S2A). We found no difference in the temporal sequence of differentiation between male and female hbNES cells, as both male (SAI1 and SAI2) and female (SAI3 and SAI5) lines exhibited similar response to directed differentiation (Fig. 2B, Fig. S2A-C). To further profile global transcriptomic changes at these differentiation timepoints, we conducted bulk RNA-seq and found that it corroborated the trends observed in our qRT-PCR results (Fig. S2B). Differentiation timepoints again tended to cluster together across cell lines, suggesting similar responses to differentiation in independently derived hbNES cell lines (Fig. S2C). Two-dimensional principal component analysis (PCA) supported this conclusion and further revealed a spectrum of progenitor states that differentiating hbNES cells transit through before becoming mature neurons (Fig. 2B). From our PCA results, we noticed one set of genes, primarily contributing to principal component (PC) 1, whose expression varied monotonically with differentiation time (Fig. 2B). However, genes contributing to PC2 had similar expression at days 0 and 56 but were maximally different at intermediate timepoints (Table S1E, Fig. S2D). We conclude from these observations that progenitor states arise and diminish during our differentiation procedure, highlighting a developmentally relevant feature of our system.

To capture the transcriptional and epigenetic changes occurring during CGN differentiation at a single cell level, we conducted scMultiome-seq at days 0, 7, 14 and 56 in a representative hbNES cell line (SAI5). This method allowed us to profile transcriptional and epigenetic landscapes from the same cell at key differentiation timepoints. We sequenced 8370, 8208, 12,952 and 7914 cells that passed quality control (Table S2A) at each timepoint, respectively. When we plotted chromatin accessibility (Fig. 2C) or integrated transcriptome (Fig. 2D) individually, we observed a continuity of cell states through days 0, 7 and 14, with day 56 situated alone. We then combined both modalities (Fig. 2E) and observed a similar trend. We noted continuity between days 0 and 7, with day 14 situated proximally and day 56 forming a distinct population. These results align with our qRT-PCR and bulk RNA-seq observations in which differentiation causes gradual changes to cell identity over time, such that day 56 cells are considerably different from those of the earlier timepoints.

To interrogate the transcriptional and epigenetic heterogeneity arising throughout differentiation, we next clustered each timepoint independently and performed differential expression (DE) and differential TF-binding motif activity (DA) analysis (Tables S2C-F, Fig. S3A,B), and computed cell cycle (i.e. S and G2/M) scores for each timepoint and cluster (Fig. S3C). Using markers reported in previously published datasets, we annotated each cluster (Table S2B). We then overlaid these clusters onto a multi-modal UMAP plot of all timepoints together (Fig. 2F). Finally, we quantified cell types at each differentiation timepoint (Fig. 2G). In our differentiation procedure, we have the emergence of an RL state at day 7, which includes human specific cell identities. This is followed by a GNP state at day 14, which matures to granule neurons at day 56 at a frequency of 86%. Interestingly, we also noted the presence of UBCs, which are thought to derive from RL in mice, at day 56 at a proportion of 5%. From these data, we conclude that our differentiation procedure is robust and developmentally relevant.

To further characterize the RL state that arose during our differentiation procedure, we first turned to our bulk RNA-seq data. We observed that known RL markers OLIG3 (Lowenstein et al., 2021), ATOH1 (Ben-Arie et al., 1997), and BARHL1 (Rachidi and Lopes, 2006) peaked in expression at day 7 (Fig. 3A-C) and were significantly upregulated at this timepoint (Table S1H). Intrigued by this finding, we turned to our scMultiome-seq data and conducted DE analysis between all timepoints. Interestingly, three of the top differentially expressed genes at day 7 relative to other timepoints were SLC12A2, TMEM132D and CRABP1 (Fig. 3D-F). SLC12A2 and TMEM132D have been described as highly expressed in the ‘Upper (Rostral) Rhombic Lip’ gene set of the Allen Brain Atlas prenatal human brain tissue gene expression dataset (Miller et al., 2014). CRABP1 expression has been identified within the chick RL (Wilson et al., 2007). Importantly, all the markers we used to characterize this timepoint as RL were either lowly or not expressed in undifferentiated hbNES cells (Fig. 3A-F). Therefore, our bulk RNA-seq and scMultiome-seq data support the emergence of an RL state at day 7 of differentiation.

To further characterize heterogeneity present within the day 7 population, we performed clustering, revealing seven subclusters, among which we performed DE analysis to identify them (Fig. 3J, Fig. S4B). Of the seven, we identified five as clusters of RL cells, all of which exhibited a high expression of SLC12A2, TMEM132D and CRABP1 (Fig. S4A). We annotated a cluster of GABAergic progenitors based on TFAP2B and PAX3 expression. We were left with a cluster of 88 cells that expressed no known granule lineage markers and proved challenging to assign an identity to. We designated this cluster Unannotated-1 (Un-1). However, based on differential expression of extracellular matrix ECM genes (DCN, COL6A3 and LUM) and the position of this cluster adjacent to hbNES cells (Fig. 2F), we suggest it represents a miniscule population of hbNES cells that are adopting a mesenchymal identity instead of differentiating.

It has recently been reported that the human RL contains a sub-ventricular compartment that is absent in mice and primates (Haldipur et al., 2019). To characterize the RL-VZ and RL-SVZ identity among our RL clusters, we adopted a gene set enrichment analysis (GSEA) strategy, using published RL, RL-VZ and RL-SVZ gene sets (Aldinger et al., 2021; Smith et al., 2022). We observed enrichment of an RL signature across all our RL clusters (Fig. 3G). Interestingly, all RL clusters showed simultaneous enrichment of RL-VZ and RL-SVZ signatures (Fig. 3H,I). A possible explanation for this observation is that in the developing human cerebellum, RL-VZ and RL-SVZ compartments are anatomically separated by a vascular plexus after 11 post-conception weeks, suggesting the requirement of such a microenvironmental cue in the specification of these distinct cell states (Haldipur et al., 2019). In the absence of such a cue, the RL may remain a mixture of these two identities. This possibility is supported by recent transcriptional profiling of human cerebellar organoids, in which RL sub-compartments were also not discerned (Nayler et al., 2021). Nonetheless, enrichment of an RL-SVZ signature within our RL population highlights the ability of our differentiation procedure to generate human-specific progenitor states.

To identify TFs that may be regulating RL identity in our system, we conducted a DA analysis of day 7 cells compared with other timepoints (Table S2L). Interestingly, two of the most active binding motifs were those of FOS (Fig. 3K) and JUN (Fig. 3L). These two TFs have recently been reported to be expressed in progenitor states of the CGN lineage, and function to prevent precocious differentiation (Goodwin et al., 2021).

Therefore, our bulk RNA-seq and scMultiome-seq data suggest emergence of an RL-like state at day 7 of differentiation at ∼92% efficiency, which includes human specific RL-SVZ identities.

RL progenitor cells migrate over the surface of the cerebellum to form a superficial layer of GNPs, which proliferate in response to SHH secreted by Purkinje neurons (Miale and Sidman, 1961; Wechsler-Reya and Scott, 1999). ATOH1 expression is maintained by GNPs until they begin to terminally differentiate (Ben-Arie et al., 1997). Given our observation of an RL-like state at day 7, we asked whether it was followed by a GNP state at day 14.

We first analyzed SHH pathway activity during differentiation. GLI1 was most highly expressed at day 14 compared with other timepoints (P=3.18e⁻¹¹) (Fig. 4E). The maintenance of ATOH1 expression at this timepoint (Fig. 3A) suggested the presence of a GNP state, which we sought to investigate further. As commercially available anti-ATOH1 antibodies are unreliable, we generated an ATOH1-EGFP knock-in allele in the SAI3 hbNES cell line using CRISPR-Cas9-mediated homology-directed repair (Ran et al., 2013) (Fig. S5). We conducted ICC using this reporter line (Fig. 4A) and observed that the proportion of EGFP⁺ cells was highest at day 14 (Fig. 4B). Although ATOH1 transcript levels peak at day 7, we posit that incomplete P2A cleavage results in ATOH1-EGFP conjugates that are protected from proteasomal degradation by elevated SHH pathway activity at day 14 (Forget et al., 2014). Importantly, the fraction of EGFP⁺ cells that were also in the cell cycle (i.e. KI67⁺) was also highest at day 14 (Fig. 4C,D). This suggests that proliferation of EGFP⁺ cells in our system is driven by SHH, which is consistent with our understanding of the effect of SHH on mouse GNPs during cerebellar development (Wechsler-Reya and Scott, 1999).

To dissect the cellular heterogeneity and identify active TFs within our GNP state, we turned to our scMultiome-seq data. We first performed DE analysis among timepoints, which further supported the GNP identity of day 14, based on upregulation of known GNP markers DCC, ZIC1 and ZIC3 (Aruga et al., 2002; Harter et al., 2010; Kaslin et al., 2009) (Table S2G). Interestingly, among the most differentially expressed genes at day 14 were ID1, ID2 and ID3: downstream effectors of BMP signaling (Ogata et al., 1993) (Table S2G, Fig. S6A). This corroborates recent data showing that BMP signaling is active in human GNPs and is necessary for EGL formation in chickens (Rook et al., 2020 preprint). We next clustered day 14 independently. We identified eight clusters, among which we performed DE analysis (Fig. 4F, Table S2E). Based on expression of known markers, we identified four clusters of GNPs (based on expression of ID1, ID2, ID3, ATOH1, DCC and DCX), two RL clusters (based on TMEM132D and CRABP1 expression) and one small glial cluster (based on SCRG1 and TSC22D expression; Aldinger et al., 2021; Lee et al., 2022) comprising 67 cells. As on day 7, we observed a small cluster of 55 cells that proved difficult to annotate that we designated Unannotated-2 (Un-2). Interestingly, Un-2 expresses the epithelial-to-mesenchymal transition (EMT) master regulator TWIST1, along with other ECM genes, suggesting mesenchymal character like the Un-1 cluster. We subsequently performed GSEA on day 14 subclusters using published RL and GNP gene sets (Aldinger et al., 2021; Smith et al., 2022). We observed enrichment of RL and GNP signatures across day 14 RL and GNP clusters, supporting our cluster identification and highlighting the biological similarities among these cell states (Fig. S6C).

At the epigenetic level, TF-binding motifs for the known GNP TFs NR4A2, MEIS1 and ZIC2 are significantly differentially active at day 14 (Fig. 4G-I), lending further support to a GNP identity being present at this timepoint (Aldinger et al., 2021; Barneda-Zahonero et al., 2012; Owa et al., 2018). Together, our bulk RNAseq, scMultiome-seq and ICC data suggest RL progenitors further mature to GNPs by day 14. GNPs comprise 84% of day 14 cells, whereas RL cells make up the remaining 16%.

As hbNES cells differentiated, we observed a gradual and concomitant decrease in the proportion of proliferating cells (Fig. S7B,C). Expression of a neuroepithelial gene signature was also gradually lost during differentiation (Fig. S8D), suggesting acquisition of a post-mitotic state and loss of stemness. Contemporaneously, differentiating cells started to extend βIII-tubulin⁺ neuronal arbors, which increased in complexity with differentiation time (Fig. S7A). Importantly, we observed few or no GFAP⁺ cells during all timepoints, suggesting that our differentiation is almost entirely neuronal.

We next sought to determine whether the cells we had obtained at our day 56 endpoint were indeed CGNs. Using our bulk RNA-seq data, we conducted DE analysis between days 0 and 56. We noticed two genes expressed by mature CGNs, SLC17A6 and GABRA6, to be significantly upregulated at day 56 (Fig. 5A,B) and we further confirmed expression of these markers at the protein level by ICC (Fig. 5C,D). We next performed GSEA of the significantly DE genes between these two timepoints identified through bulk RNA-seq (Fig. S8A, Table S1M-O). We noticed several pathways pertaining to glutamatergic synaptic transmission were significantly enriched at day 56 relative to day 0 (Fig. S8A). Notably, one of these terms was ‘parallel fiber to Purkinje cell synapse’, which refers to the glutamatergic synapse formed by CGN axons onto Purkinje neuron dendrites in the cerebellar molecular layer (Fig. S8B). These results suggest that the cells we obtain at the end of our differentiation have a mature CGN identity.

We further characterized the day 56 population at the single cell level. Relative to other timepoints, CGN markers RELN, GRIK2 and NFIA (Aldinger et al., 2021; Goldowitz et al., 1997) were significantly upregulated at day 56 (Table S2G). Enrichment of a granule neuron signature was also highest at day 56 (Fig. S8D). We then clustered day 56 independently and discovered six clusters, among which we performed differential expression analysis (Table S2F). We were able to annotate these clusters based on the expression of known markers as CGNs (based on expression of ROBO2, RELN, GRM8, NFIB and FSTL4; Aldinger et al., 2021), UBCs (TRPM3 and OTX2) (Vladoiu et al., 2019) or RLs (TOP2A, ASPM and CENPE) (Aldinger et al., 2021) (Fig. 5E, Table S2B). The overall output of our differentiation procedure was 86% granule neurons, 8% RL and 6% UBCs. We corroborated the presence of UBCs at this timepoint by conducting differential TF motif activity analysis among day 56 clusters. In addition to OTX2 expression (Fig. S9A), we found significantly enriched motif activity for OTX2 and EOMES in the UBC cluster (Fig. S9B,C). As UBCs are thought to arise from RL progenitor cells (Englund et al., 2006; Vladoiu et al., 2019), we were curious about whether this was the case for our system. We applied transcriptional lineage reconstruction using URD (Farrell et al., 2018), and observed that a likely bipotent RL cell generates the UBCs we observe at day 56 (Fig. S9D). Although UBCs have been previously generated from human pluripotent stem cell-derived cerebellar organoids (Nayler et al., 2021; Silva et al., 2021), we extend these observations by suggesting a developmental origin for these cells. The persistence of an RL population at day 56 confirms latest reports of human cerebellar development, in which the RL is a long-lived germinal zone, present for up to the first 2 years of life (Haldipur et al., 2019); we further confirmed the identity of this day 56 RL population by GSEA (Fig. S8C).

To assess the functional potential of our model, we performed electrophysiological measurements of neuronal activity. We differentiated hbNES cells on microelectrode arrays (MEA) to record spontaneous extracellular activity at a population level every week from days 14 to 56. Extracellular spikes were first detected on day 28 and continued to be recorded with increasing frequency until day 56 (Fig. 5F, Fig. S7D), representing functional neuronal activity. To determine whether glutamate was the neurotransmitter responsible for driving the spontaneous activity of our neurons, we pharmacologically blocked all post-synaptic glutamate receptors with the small molecules CNQX (20 µM) and D-AP5 (100 µM), causing a significant decrease in mean firing rate (Fig. 5G). We next assessed neuronal function at the single cell level using current-clamp whole-cell electrophysiology. We observed that neurons at day 56 can generate single evoked action potentials (Fig. 5H); however, they remain depolarized after stimulation and do not fire repeatedly. This can be attributed to a combination of lack of supporting cell types, absence of neural circuit integration (Mattugini et al., 2019) and absence of appropriate voltage-gated potassium ion channels required to sufficiently hyperpolarize the cell membrane. We quantified recovery from after-hyperpolarization (Tau_AHP) on neurons we patched (Fig. 5I) and observed that they were in agreement with Tau_AHP values previously described for human ES cell-derived neurons (Schrenk-Siemens et al., 2015). Taken together, our data indicate that the predominant output of our differentiation protocol is CGNs that are capable of spontaneous and evoked neural activity in vitro.

Having generated and characterized a system of CGN development, we were next interested in studying the role of SHH in this process. Previous reports have shown that SHH acts as a mitogen for GNPs (Wechsler-Reya and Scott, 1999) and is also required for cerebellar GABAergic neural development in mice (Huang et al., 2010). Therefore, we hypothesized that SHH plays a similar role in human CGN differentiation. We designed an experimental approach in which we exposed differentiating hbNES cells to three different concentrations of SHH in the growth medium starting at day 7 (Fig. 6A). We noticed a dose-responsive relationship between SHH target GLI1 expression and SHH concentration, verifying our experimental approach (Fig. 6B). Interestingly, we observed that SHH does not appear to be required for commitment to the granule lineage, as ATOH1 expression was unchanged among all three SHH conditions (Fig. 6C). This result reinforces the idea that granule lineage commitment likely occurs at an RL stage in the human context before progenitor cells migrate to form the EGL. Furthermore, we noted a dose-responsive relationship between SOX2 expression and SHH concentration (Fig. 6D). Previous work demonstrating that GLI transcription factors regulate SOX2 expression provides a possible mechanism for this observation (Bora-Singhal et al., 2015), and is consistent with our previous identification of a SOX2⁺ population persisting in the mouse EGL (Selvadurai et al., 2020). Next, we sought to investigate the role of SHH as a mitogen during human CGN differentiation. Using our ATOH1-EGFP reporter line, we performed KI67 staining at day 14 of differentiation to quantify the fraction of proliferating cells. We observed that the fraction of GFP⁺ cells was similar under all three concentrations of SHH (Fig. S10A,C), supporting our qRT-PCR result. On day 14, we observed a positive correlation between the concentration of SHH in the growth medium and the proportion of GFP⁺ cells that were also KI67⁺ (Fig. 6G,H). This effect was also observed at the population level across the timepoints we analyzed (Fig. S10B). This result further confirms the physiological relevance of our system, as it demonstrates that SHH acts as a mitogen of GNPs.

To test whether SHH promotes GABAergic fates in the CGN lineage, we analyzed the expression of SLC17A7 and GAD1, glutamatergic and GABAergic markers, respectively, during differentiation using qRT-PCR (Fig. 6E,F). Neurons derived under all three conditions showed similar extents of arborization (Fig. S10E). However, we observed that the GABAergic marker GAD1 was significantly upregulated in neurons differentiated in 15 ng/ml of SHH compared with 5 ng/ml (Fig. 6F). SLC17A7 expression was also decreased in neurons differentiated under high SHH, although not significantly (Fig. 6E). We then asked whether these fate differences manifest physiologically. Using MEA methodology, we observed that although functional neurons were generated under all conditions (Fig. 6I, Fig. S10F), those exposed to moderate levels of SHH had the highest levels of spontaneous activity (Fig. 6J). Therefore, we propose a model whereby high SHH activity potentiates an immature cell state and the GABAergic lineage, leading to less spontaneous neuronal activity at the population level. However, in the absence of SHH, functional neuronal maturation is also impaired, perhaps suggesting a pro-differentiation role of SHH within a specific concentration range.

Diseases of the granule neuron lineage arise in progenitor states for which human models have hitherto been lacking. Through the directed differentiation of human hbNES cells, we present a method to efficiently obtain RL and GNP cell states. Extended differentiation of hbNES cells produces functional mature CGNs at 56 days. Importantly, we observed recapitulation of a human specific RL-SVZ compartment, and persistence of an RL identity to day 56. These findings are consistent with recent reports of maintenance of an RL progenitor pool for as long as the first 2 years of life (Aldinger et al., 2021; Haldipur et al., 2019).

Although the CGN lineage was the predominant output of our differentiation, we observed small populations of GABAergic progenitors and glia at days 7 and 14, respectively. The likeliest explanation for these populations is that a minority of very strongly GABAergically primed hbNES cells commits to the cerebellar ventricular zone lineage upon EGF and FGF withdrawal. We propose that such cells first generate GABAergic VZ progenitors at day 7 and then glia at day 14, which are known to derive from the VZ. We also observed small clusters of cells at days 7 and 14, which proved difficult to annotate but expressed markers suggestive of mesenchymal identity. Given the positions of these clusters adjacent to undifferentiated hbNES cells on our global UMAP (Fig. 2F), we speculate that they may represent small populations of hbNES cells that were resistant to differentiation and instead adopted mesenchymal character.

Despite the absence of a three-dimensional environment and lack of interaction with other supporting cell types that facilitate CGN differentiation in vivo (Consalez et al., 2021; Leto et al., 2016), we nonetheless recapitulate the full temporal sequence of human CGN development from hbNES cell to functional neuron. By providing an in-depth resource of human CGN development at the single-cell resolution, our work complements recent studies of human cerebellar development using organoids and primary tissue (Aldinger et al., 2021; Nayler et al., 2021; Silva et al., 2021). We present a new method by which hypotheses pertaining to CGN development and disease can be functionally tested.

The ability to isolate and study the contributions of individual growth factors is important to understanding developmental processes. Such experiments are technically challenging to perform in vivo, because it may not be feasible to specifically modulate the level of growth factor a certain cell type is exposed to in a whole organism. Additionally, individual growth factors may regulate several developmental processes simultaneously, which may confound the phenotype of interest. Therefore, in vitro systems such as ours are best suited to answer such questions. We leveraged this strength of our system to test the role of SHH during CGN differentiation, revealing an intriguing role of SHH in lineage determination and functional maturation. Our observation that elevated SHH signaling potentiates an immature cell state and inhibits functional neuronal maturation is consistent with the finding that constitutive SHH pathway activation can transform granule lineage cells to seed medulloblastoma (Schüller et al., 2008; Selvadurai et al., 2020; Vanner et al., 2014; Yang et al., 2008).

The ability of embryonically derived hbNES cells to stably maintain their regional and temporal identities in vitro represents perhaps the most powerful aspect of this system. As all cerebellar tissue arises from r1, it is conceivable that with application of the appropriate combination of growth factors for an optimal amount of time, all human cerebellar lineages can be derived in vitro, thereby generating experimental models for the full spectrum of cerebellar pathologies.

An exciting observation of our day 56 population was the emergence of UBC-like cells. UBCs are glutamatergic interneurons found in the cerebellar granule layer. To our knowledge, our work is the first evidence that human UBCs derive from RL progenitor cells, as they do in mice. Our work thereby lays the groundwork for further studying the development of this cell type in the human context. Considering recent reports implicating UBC lineage stem and progenitor cells as origins for group 3 and group 4 medulloblastoma (Hendrikse et al., 2022; Smith et al., 2022), further study of this lineage and its use in tumor modeling will be of great clinical interest.

Overall, we define the complete temporal sequence of human CGN lineage development in vitro, highlighting the cell state transitions from a neuroepithelial stem cell population to committed neural progenitors, to differentiated and functional granule neurons. By defining the distinct steps of lineage transition, this system provides an opportunity to model the timing and cellular context of disease-associated alterations and their effects on human cerebellar cellular maturation.

All experiments with human cells were approved by The Hospital for Sick Children Research Ethics Board (REB 0020020272). hbNES cells isolated from primary human fetal sources at gestational weeks 5 to 6 (Carnegie stages 15-17) (Tailor et al., 2013) were grown in DMEM/F12 (Wisent, 319-075-CL), 1× N2 (homemade), 0.01× B27 (Life Technologies, 17204-44), 10 ng/ml rhEGF (Sigma, E9644) and 10 ng/ml bFGF (STEMCELL Technologies, 02634) in a humidified 37°C, 5% CO₂ incubator on Corning Primeria or Falcon tissue culture vessels coated with poly-L-ornithine (Sigma, P4957-50ML) and laminin (Sigma, L2020-1MG). hbNES cell lines tested negative for Mycoplasma. When cells neared confluency, they were dissociated with TrypLE express (Life Technologies, LS12605028) for 5 min. The dissociation was quenched with DMEM/F12, cells were pelleted by centrifugation at 400 g for 5 min and then replated at desired density. Cells were cryopreserved in growth medium supplemented with 10% DMSO at −80°C or in liquid nitrogen. Growth medium was partially replaced every 24 h.

A guideRNA sequence directing Cas9 cleavage 2 bp away from the desired site of knock-in was found using Benchling (5′-TGACTCGGATGAGGCAAGTT-3′; on target score, 25.8; off target score, 75.3). Sense (CACCGTGACTCGGATGAGGCAAGTT) and anti-sense (AAACAACTTGCCTCATCCGAGTCAC) oligonucleotides were hybridized and ligated into pSpCas9(BB)-2A-GFP (PX458) (Addgene plasmid 48138; a gift from Feng Zhang) as previously described (Ran et al., 2013). A homology repair template was designed and ordered from VectorBuilder. hbNES cells were nucleofected using the Amaxa 4D-Nucleofector X-unit system (Lonza Bioscience, AAF-1002X) when they had grown to 75% confluency. 2×10⁶ cells were pelleted and resuspended in 100 µl of P3 solution containing 2.5 µg of PX458 and 2.5 µg of homology template. The cell suspension was transferred into a 100 µl nucleofection cuvette (Lonza Bioscience, V4XP-3024) and nucleofected with code DN100, which has previously been demonstrated to deliver plasmid cargo with high efficiency to adherent neural stem cells (Bressan et al., 2017). After nucleofection, cells were transferred to one well of a six-well plate and allowed to recover for 6 days. Cells that had undergone homologous recombination were selected with 1 mg/ml G418 (Invitrogen, 11811023). The presence of the knock-in allele was detected by PCR (see Table S1B for primers). Homologous recombinants were further isolated by twice sorting the population for mCherry⁺/DAPI⁻ cells on a Sony MA900 cell sorter. Genotyping PCR primers are listed in Table S1B.

6×10⁵ hbNES cells were plated in growth medium on poly-L-ornithine/laminin-coated 6 cm plates on the day before differentiation, and 2.5×10⁴ hbNES cells were plated on coated glass coverslips. On day 0 of differentiation, growth medium was completely aspirated and replaced with Step 1 differentiation medium: DMEM/F12, 1× N2, 0.01× B27, 2.5 ng/ml BMP6 (R&D Systems, 507-BP-020/CF), 5 ng/ml BMP7 (R&D Systems, 354-BP-500/CF) and 50 ng/ml BMP12/GDF7 (R&D Systems, 8386-G7-050). On day 4, differentiating cells were fully confluent and split 1:3 into Step 1 differentiation medium. On day 7, differentiating cells were confluent again and split 1:3-1:4 into Step 2 differentiation medium: BrainPhys Neuronal Medium (STEMCELL Technologies, 05790), 1× N2, 0.01× B27, 2.5 ng/ml BMP6, 5 ng/ml BMP7, 50 ng/ml BMP12/GDF7, 5 ng/ml SHH (R&D Systems, 8908-SH), 2 ng/ml JAG1 (1277-JG-050, 1277-JG-050), 10 ng/ml NT3 (R&D Systems, 267-N3-025/CF) and 100 ng/ml BDNF (R&D Systems, 248-BDB-01 M/CF). From days 0 to 28, medium was partially replaced every 48 h. After day 28, medium was partially replaced twice per week. For the SHH dependence experiment, two additional formulations of step 2 differentiation medium were prepared with 15 ng/ml SHH or 0 ng/ml SHH.

Total RNA was extracted from cell pellets using the Qiagen RNeasy Plus Mini Kit (Qiagen, 74134). cDNA was synthesized with Transcriptor Reverse Transcriptase (Roche, 3531287001), using an oligo(dT)₁₈sequence. Reactions were carried out in either 96-well (Bio-rad, HSP9661) or 384-well (Thermo Fisher Scientific, 43-098-49) format using SsoFast EvaGreen Supermix (Bio-rad, 1725202) or PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, A25778), respectively. Primer specificity and amplification efficiency for genes of interest were verified by performing reactions on a dilution series of human reference cDNA (Takara, 636693). Gene expression relative to GAPDH and RPLP0 house-keeping genes was calculated using Microsoft Excel, and data were visualized using GraphPad Prism and Morpheus.

Cells grown on glass coverslips were fixed in 10% formalin for 10 min at room temperature. Cells were blocked and permeabilized in a PBS solution of 0.1% Triton X-100 and 5% normal goat serum (NGS) for 1 h. Cells were incubated with primary antibody at the indicated concentration for either 1 h at room temperature or overnight at 4°C, washed three times, and then incubated with secondary antibody for 1 h at room temperature. See Table S3 for further information on primary and secondary antibodies used. Coverslips were mounted on glass slides with fluorescence mounting medium. (Dako, 34538). Cells were imaged on a Zeiss Epifluorescence or Leica SP8 Lightning STED confocal microscope at the SickKids Imaging Facility. Images were processed in Fiji and cells were counted using the built-in cell counter plug-in (Schindelin et al., 2012). At least three biological replicates of each staining were counted, and analyses were unified by cell count.

24- or 48-well CytoView plates with 16 recording channels per well were obtained from Axion Biosystems (M384-tMEA-24W and M768-tMEA-48W). MEA plates were coated with poly-L-ornithine overnight, washed four times with water, and then coated with laminin for at least 24 h before cell plating. On day 7 of differentiation, 1.5×10⁴ cells were plated into each well in a 70 µl droplet, containing laminin at a 1:100 dilution. Cells were allowed 60 min to adhere after which medium was topped up to 200 µl (48-well plate) or 500 µl (24-well plate). MEA plates were recorded for spontaneous activity every 7 days from day 14 of differentiation to day 56 using the AxIS software (version 2.5.2.1) on the Maestro Original system (Axion Biosystems). Recordings were performed at a constant temperature of 37°C, and plates were allowed to rest undisturbed in the pre-warmed reader before spontaneous neural activity was recorded for an additional 5 min. Pharmacology experiments were performed at day 56. CNQX (Biotechne, 0190) and D-AP5 (Biotechne, 0106) were reconstituted in DMSO and artificial cerebrospinal fluid (ACSF), respectively. Cells were first equilibrated in growth medium containing vehicle (i.e. DMSO and D-AP5) for 10 min before activity was recorded for 5 min. Medium was then completely removed and cells were washed once with PBS before being equilibrated in growth medium containing CNQX (20 µM) and D-AP5 (100 µM) for 10 min, before activity was recorded for 5 min. Recording channels were sampled at a frequency of 12.5 kHz and bandpass filtered from 200 Hz to 3 kHz. Spikes were then detected using an adaptive threshold-crossing method, with the threshold set to six times the standard deviation of the estimated noise of each channel. Offline analysis was conducted using the Neural Metric Tool (Axion Biosystems, version 3.1.7). Weighted mean firing rate (wMFR) was calculated as the mean firing rate normalized by the number of active electrodes in each well. Electrodes were considered active if at least five spikes per minute were detected. Single channel burst detection was performed using the Poisson surprise algorithm, with the ‘minimum surprise’ parameter set to 3. Spike and burst metrics obtained using the Neural Metric Tool were subsequently visualized using GraphPad Prism. Extracellular action potential waveforms were obtained from AxIS software and visualized in GraphPad Prism.

NES cells differentiated to day 56 on glass coverslips were transferred to a recording chamber filled with bath solution. The bath solution consisted of (in mM) 140 NaCl, 5 KCl, 2 CaCl₂, 2 MgCl₂, 10 glucose and 10 HEPES (pH adjusted to 7.4 with NaOH). Patch pipettes for recording, with a resistance of 3-4 MΩ, were filled with intracellular solution consisting of (in mM) 123 K-Gluconate, 10 KCl, 1 MgCl₂, 1 EGTA, 10 HEPES, 0.1 CaCl₂, 1.2 ATP and 4 glucose (pH adjusted to 7.2 with KOH). Whole-cell recordings were made with a MultiClamp 700B amplifier paired with a Digidata 1440A digitizer and acquired online using the pClamp10 software package. In current-clamp experiments, cells were injected with minimal current to maintain their voltage at around −60 mV and current steps were injected starting with −50 pA for 500 ms, increasing each sweep by 50 pA to a final amount of 650 pA. Representative trace in the figure demonstrates that the action potential (AP) was elicited at 50 pA current injection. All experiments were performed at room temperature. AP traces were quantified and graphed using GraphPad Prism and Adobe Illustrator.

250 ng of RNA from each sample were used for library preparation using the NEBNext Ultra II Directional poly(A) mRNA library prep kit (E7760S). Paired-end sequencing (2×100 bp) was then conducted on an Illumina Novaseq S1 flowcell. Reads were mapped to release 36 (GRCh38.p13) of a human reference transcriptome using Salmon v1.1.0 (Patro et al., 2017). Downstream analysis was performed in R using tximport (Soneson et al., 2016) and transcript abundances were log-normalized before differential expression analysis with the DESeq2 package (Love et al., 2014). For DE analysis, a log fold change of 1 and a false discovery rate of 0.05 were set. P-values were calculated using the Wald test and were adjusted for multiple hypothesis testing. Gene set enrichment analysis was performed using the ClusterProfiler package in R (Wu et al., 2021).

Cryopreserved cell suspensions were prepared in growth medium supplemented with 10% DMSO. Cells were thawed and nuclei were isolated according to protocol CG000365, Revision C by 10X Genomics (Nuclei Isolation for Single Cell Multiome ATAC+Gene Expression Sequencing). After nuclear isolation, GEMs were generated and barcoded using the Chromium Next GEM Chip J. ATAC and gene expression libraries were then constructed according to protocol CG000338, Revision E by 10X Genomics (Chromium Next GEM Single Cell Multiome ATAC+Gene Expression). Library size was subsequently quantified using a BioAnalyzer. All samples were sequenced on the Illumina NovaSeq 6000 system, with a target of 5×10⁴ RNA reads/nucleus and 2.5×10⁴ ATAC reads/nucleus. Sequencing reads were aligned to the human reference genome and pre-processed using the CellRanger v2.0.0 pipeline.

Downstream analysis of scMultiome data were conducted using the Seurat and Signac packages in R (Hao et al., 2021; Stuart et al., 2021). Low-quality cells were first filtered out according to metrics reported in Table S2A. Doublets were identified using the DoubletFinder package with the following parameters: pK=0.1, pN=0.05 (day 0), 0.01 (day 7), 0.005 (days 14 and 56), nExp_poi=0.075 * number of cells in object, and then removed (McGinnis et al., 2019). Next, peaks were called using MACS2. Gene expression data for individual timepoints were then normalized using SCTransform. All four timepoints were integrated in Seurat across their ‘RNA’ assays, before gene expression was scaled. Dimensionality was reduced with principal component analysis (PCA). DNA accessibility data were processed by latent semantic indexing (LSI). A multimodal nearest neighbor graph was generating using the top 30 PCs and top 40 LSI, after which a joint UMAP plot was constructed to visualize both data modalities. Individuals timepoints were clustered at a resolution of 0.2, and DE analysis was conducted between clusters using a log2 fold change of 0 and a maximum of 1000 cells per identity using a Wilcoxon Rank Sum test; multiple hypothesis testing was corrected for. DE was also conducted between timepoints using the same parameters. Cell clusters were annotated using markers identified through DE (see Table S2B). Differential accessibility (DA) analysis was performed between day 0 and days 7, 14 and 56 using the above parameters. DA regions in each differentiation timepoint relative to day 0 with P<0.005 were used to search for transcription factor motifs present in these regions. Motif activity scores were then computed using chromVAR (Schep et al., 2017).

The escape package in R was used to performed single cell GSEA (Borcherding et al., 2021). RL-VZ, RL-SVZ, EGL and UBC gene sets were obtained from Smith et al. (2022) and RL, GNP and GN gene sets were obtained from Aldinger et al. (2021). The NES gene set consisted of: LIN28A, LIN28B, PRTG, SOX1, SOX2, ADGRL2, ADGRL4 and ROR2.

Transcriptional reconstruction of the UBC lineage was performed using URD (Farrell et al., 2018), and the transcriptomes of cells from days 7, 14 and 56. RL-1 was designated as the root cluster and GN-1 and UBC were set as the tip clusters.

No custom software was written for data analysis. Raw and processed sequencing data have been deposited in GEO under accession number GSE233339.

Data were analyzed using GraphPad Prism 9. Two-tailed parametric t-tests were performed to compute P-values, and a threshold of P<0.05 was used for statistical significance. Data are mean±s.e.m., unless otherwise specified.

We thank Prof. Austin Smith for the hbNES cells used in this study. We acknowledge support from The Centre for Applied Genomics (TCAG) at The Hospital for Sick Children for their bulk RNA-seq service and Princess Margaret Genomics Centre (PMGC) for their scMultiome-seq service.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201534.reviewer-comments.pdf