Single-cell transcriptomics of the early developing mouse cerebral cortex disentangle the spatial and temporal components of neuronal fate acquisition

In the nervous system of vertebrates, the production of assorted types of neurons enables the execution of elaborate cognitive functions. In recent years, single-cell RNAseq (scRNAseq) approaches revealed that neuronal diversity is far more extensive than anticipated (Tasic et al., 2018). The developmental mechanisms enabling the emergence of such a diversity are currently under intense investigation.

The neocortex, which has been subjected to strong expansion of cell numbers in primates and especially humans, has attracted most of the attention so far. It is organised in six layers of neurons that are generated according to an ‘inside-out’ sequence (Angevine and Sidman, 1961) and is flanked by three layers of cortices: the hippocampus medially and the piriform cortex laterally. All glutamatergic cortical neurons are derived from apical progenitors (APs) that reside in the pallial ventricular zone (VZ) and divide to self-renew, to produce neurons (direct neurogenesis) or to give rise to basal progenitors (BPs), which display a limited proliferation potential in mice (indirect neurogenesis). In the subpallium, ventricular progenitors give rise to inhibitory neurons that will either populate the striatum or migrate dorsally to invade the cortical plate. The so-called pallial-subpallial boundary (PSB) appears initially as a fuzzy border between the dorsal and ventral aspects of the developing telencephalon that is progressively refined to precisely separate the two compartments (Cocas et al., 2011; Corbin, 2003). The PSB not only delimitates distinct regions but also acts as a secondary signalling centre, releasing morphogens crucial for cortical patterning (O'Leary and Sahara, 2008). Within the pallium itself, no sharp boundaries between progenitor domains can be found. Instead, the graded expression of patterning genes, such as Pax6, Lhx2 and Emx1/2, led to the proposal that diversity is not strictly spatially segregated at the progenitor level, as is the case in the developing spinal cord (Pierani and Wassef, 2009; Sagner and Briscoe, 2019; Yun et al., 2001). How the apparent continuity of progenitors gives rise to discrete neuronal types has not been clearly established. Currently, temporal regulations are believed to contribute most to neuronal diversity in the cerebral cortex. According to this model, the neurogenic differentiation followed by neocortical neurons is essentially the same and it is the temporal maturation of APs that will influence the final neuronal output (Okamoto et al., 2016; Telley et al., 2019). However, temporal regulations are unlikely to account for all the neuronal diversity in the developing pallium, otherwise the hippocampus and piriform cortex would not significantly differ from the neocortex. The tetrapartite model of pallial development (Puelles et al., 1999, 2000) postulates that APs segregate in four consecutive domains that are conserved among vertebrates: the medial pallium, giving rise to the hippocampus; the dorsal pallium (DP), giving rise to the neocortex; and the lateral and ventral pallium (LP and VP), contributing to the piriform cortex, claustrum and amygdala. However, the precise contribution of LP and VP domains to mature telencephalic structures are still being discussed (Watson and Puelles, 2017; Wullimann, 2017), in part due to complex migratory patterns and to the existence of both laminar and nuclear derivatives. Furthermore, the diversity and identity of LP- and VP-derived neurons remains poorly characterised.

To determine the extent of diversity among progenitors and neurons, and apprehend how discrete neuronal types arise from continuous progenitor domains, we implemented a scRNAseq approach. We first characterised neuronal diversity around the PSB, mapped the identified cell types on the tissue and assessed their ontogeny using genetic tracing. We then investigated AP diversity and extensively characterised the molecular signatures of progenitor domains on both sides of the PSB. We used previously published data exploring the temporal maturation of dorsal cortex progenitors (Telley et al., 2019) to discriminate genes whose expression in progenitors changes with the developmental stage from those subjected to spatial regulations along the dorso-ventral (DV) axis. Finally, we reconstructed the developmental trajectories leading to the production of VP- and DP-derived neurons to identify gene waves specific to each lineage. Taken together, our result allows us, for the first time, to unravel the respective contribution of temporal maturation and spatial information to the establishment of neuronal diversity in the developing dorsal telencephalon.

Using a droplet-based approach (10x Genomics Chromium V2), we performed unbiased scRNAseq on dissected explants containing the PSB and adjacent tissue throughout the entire rostro-caudal (RC) extent of the telencephalon, obtained from wild-type E12.5 embryos (Fig. 1A). After quality control and filtering, we obtained 4225 cells sequenced at a median depth of 11,665 UMI and a median detection of 3571 genes per cell.

Visualisation of transcriptional heterogeneity on a SPRING plot revealed that the dataset is organised in a ‘crab-shaped’ topology consisting in a large group of cycling progenitors [Sox2 and Nes (Nestin)] from which stemmed two branches of postmitotic cells corresponding to pallial glutamatergic [Tbr1 and Slc17a6 (Vglut2)] and subpallial GABAergic [Dlx5 and Gad2 (GAD65)] differentiation trajectories, respectively (Fig. 1B,C). Gene waves specific or common to these lineages are detailed in Fig. S1A and Tables S2-S4.

To ensure that we successfully captured all states through which pallial cells transit from progenitor to mature neuron, we calculated expression scores for apical (AP) or intermediate/basal (BP) progenitors and neurons at early (EN) or late (LN) stages of differentiation using sets of genes previously shown to characterise these distinct steps (see Materials and Methods). Projection of the scores on the SPRING plot recapitulates the known sequence of transcriptional states by which cells transit from AP to LN, indicating that we sampled the different intermediates along the differentiation trajectory of pallial neurons (Fig. 1D).

When focusing on cells from the glutamatergic branch with high BP signature score (Fig. S1B), we noticed that one group expressing genes associated with proliferation, such as Top2a, Mki67, Cdk1 or Aurka, and annotated as S- or G2/M-phase cells, appears to diverge from another group of aligned G0/G1 cells (Fig. 1B, Fig. S1B). The latter most likely corresponds to the trajectory followed by cells transitioning from AP to neurons without achieving a last division as BP; hence, representing a non-proliferative intermediate state. This observation suggests that, at E12, direct neurogenesis involves a transition through a ‘BP-like’ transcriptomic state Eomes⁺ (Tbr2), Btg2⁺ (Tis21) and Neurog2⁺ without replicating their DNA or dividing. Taken together, these results indicate that the two main axes in the low-dimensional manifold of our dataset reflect (1) the differentiation status, from cycling progenitor to postmitotic mature neuron; and (2) the pallial or subpallial origin of the cells (Fig. 1E), recapitulating global tissue organisation around the PSB at E12.5.

To explore neuronal diversity, we first selected mature neurons based on their high LN and low EN scores, and performed iterative graph-based clustering using the scrattch.hicat pipeline (Tasic et al., 2016) to identify glutamatergic and GABAergic cell types. Three main GABAergic subtypes were identified and could be further subdivided (Fig. S2). A first group displayed the characteristic Lhx6/Sst/Calb1 signature of cortical interneurons migrating from the medial ganglionic eminence and preoptic area. Two clusters contained neurons expressing the dorsal lateral ganglionic eminence marker gene Sp8 (Ma et al., 2012), and could be split by that of either Six3 or Zic1 (Fig. S2). The remaining clusters were characterised by the strong expression of Isl1, Ebf1 and Zfp503 (Nolz1), genes that are reported to be expressed by striatal projection neurons (Chang et al., 2004; Corbin, 2003; Garel et al., 1999).

To accurately characterise glutamatergic neurons, we retained only core cells to define clusters (Fig. 2A-C). Analysis of the higher hierarchical relationship revealed three major groups (Fig. 2D) corresponding to Cajal-Retzius (CR) cells, the only Foxg1-negative cells, branching apart from two other groups distinguished by the expression of Neurod6 or Pbx3. The genes found significantly enriched in each cluster are listed in Table S5.

Consistent with recent reports, CRs were found to form a transcriptionally distinct group of neurons in the cortex (Tasic et al., 2018; Yao et al., 2021). They were readily identified by the absence of Foxg1 (Hanashima, 2004) but also Mn1, and the expression of previously described markers such as Reln, Lhx5, Calb2 (calretinin), Dmrta2 (Dmrt5), Nr2f2 (COUP-TFII) and Ebf2 (Chuang et al., 2011; Miquelajauregui et al., 2010; Ogawa et al., 1995; Saulnier et al., 2013; Tripodi, 2004) (Fig. 2D, Fig. S3B). Of note, none of these genes is fully specific to CRs; most noticeably, Reln was found in both Bhlhe22⁺ and Foxp2⁺ populations, although at lower levels than in CRs (Fig. 2D).

CRs divide in two main groups: the first one expressing a set of genes among which Trp73, Ebf3, Cdkn1a, Lhx1 and Cacna2d2; and a smaller set displaying higher, albeit not restricted, expression of Unc5c and Grm1 (encoding mGluR1; recognised by lot1 antibody), and high levels of Calb2 and Meis2 (Fig. 2D, Fig. S3I). CRs originate from several distinct regions of the forebrain and Trp73 expression is a well-described feature of those derived from medial sources, i.e. the cortical hem, pallial septum and thalamic eminence (Griveau et al., 2010; Meyer et al., 2002; Ruiz-Reig et al., 2017; Tissir et al., 2009), whereas VP-derived CRs were previously shown to be p73 negative and Calretinin^high (Griveau et al., 2010; Hanashima et al., 2007). We therefore concluded that the first hierarchical distinction between CRs subtypes is their medial versus VP origin.

Among Trp73⁺ CRs, a large cluster displayed a strong signature consisting of known genes expressed in CRs, such as Grm2 (mGluR2) or Cxcr4 (Stumm et al., 2003; Yamazaki et al., 2004), combined with dorsal cortex markers (Nfib and Nfix) (Plachez et al., 2008), whereas a smaller cluster had a Insm2⁺/Pbx3⁺/Mab21l1⁺ signature and shared the expression of Chl1 and Grm1 with VP-derived CRs (Fig. 2D, Fig. S3J).

In situ hybridisation patterns (Fig. 3A,C) indicated that Trp73-negative CRs (VP derived) are mostly confined to the lateral cortex and lateral olfactory tract (LOT) region, consistent with their low Nfib expression and with previous reports (Bielle et al., 2005). Regarding Trp73⁺ CRs (those from medial sources), the Insm2⁺/Pbx3⁺ cells were found around the LOT and in the subpallial marginal zone (MZ), whereas Grm2⁺/Nfib⁺ subtypes were found above the entire dorsal cortex (Fig. 3A,C).

We performed fate-mapping experiments in order to assess whether the two Trp73⁺ CR subpopulations correspond to distinct origins. Using the Wnt3a^Cre line to label hem-derived CRs, we found all of them to be Nfib positive, including those that migrated down close to the LOT region (Fig. 4A). By contrast, the ΔNp73^Cre line, which labels CRs from the hem, septum and thalamic eminence, allowed the identification of some Nfib-negative CRs at all medio-lateral levels (Fig. 4B). We also found Meis2-positive cells, especially around the LOT and occasionally in the DP (Fig. 4C). We therefore concluded that the Trp73⁺/Nfib⁺ subset corresponds to hem-derived CRs, whereas the Trp73⁺/Pbx3⁺ cells most likely derive from the septum and/or thalamic eminence. Our data are therefore consistent with previous histological characterisation of CR subtypes and demonstrate that VP-derived CRs clearly diverge from other subtypes generated at medial sources.

Among the Neurod2⁺/Neurod6⁺/Tcf4⁺ clusters (node G in Fig. S3A), the first population to branch apart was characterised by the specific expression of Fezf1 (Fig. 2D). In situ hybridisation for Fezf1 stained a sharp domain in the ventral aspect of the prospective insular cortex (Fig. 3B,D). The remaining two clusters shared the expression of Nfia, Nfib, Pcp4 and Fezf2, and were distinguished by expression of either Nfix and Bhlhe22 (Bhlhb5) or Ppp1r14c and Scg2 (Fig. 2D, Fig. S3H). In situ hybridisation revealed that Bhlhe22 and Nfix are expressed in the preplate of the DP, whereas Ppp1r14c and Scg2 are found more laterally (Fig. 3B,D, Fig. S4). The three Neurod6⁺ populations therefore occupy distinct dorsoventral positions. We then used genetic tracing in the Emx1^Cre line, which recombines mostly in LP and DP progenitors (Gorski et al., 2002), to determine that 96±1% of Nfib⁺ and 93±1% of Bhlhb5⁺ cells derive from Emx1 lineages and are therefore LP/DP derivatives, consistent with their position (Fig. 4D). We also found that the in situ hybridisation signal for Fezf1 overlaps with the ventral-most GFP⁺ area in Emx1^Cre;Rosa26^YFP embryos, suggesting that this population also derives from the LP/DP (Fig. 4E). However, in the absence of a LP- or DP-specific Cre driver, and despite clear distinctions in the relative DV position of Bhlhe22⁺, Ppp1r14c⁺ and Fezf1⁺ populations, it is difficult to ascertain whether they can be classified as DP or LP derivatives.

Four populations were found among the Pbx3⁺ and Mab21l1⁺ branch (Fig. 2D). Three of them showed expression of the transcription factor Foxp2 and could be distinguished by the expression of either Lhx5, Mef2c or Etv1 (Er81), although these genes were also detected in other clusters (Figs 2D and 3B,D, Fig. S3). Tissue mapping indicated that Etv1⁺ cells are found in the most rostral regions (Fig. S4, Zimmer et al., 2010). By contrast, the Lhx5⁺ population is distributed more caudally and could correspond to neurons migrating from the thalamic eminence towards the olfactory bulb (Huilgol et al., 2013). The Foxp2/Mef2c population could not be precisely mapped due to the absence of a specific marker. The last cluster corresponded to neurons specifically expressing Nr4a2 (encoding Nurr1) and Cck (Fig. 2D). Among pallial derivatives, these cells occupy the most ventral position (Fig. 3B,D) and also showed the lowest Tbr1 but highest Lhx9 expression (Fig. 2D), which is the signature of the previously described VP migrating stream (Tole et al., 2005). To determine the origin of Nr4a2⁺ and Foxp2⁺ neurons, we took advantage of the specific expression of the transcription factor Dbx1 in VP progenitors (Bielle et al., 2005; Medina et al., 2004) to perform lineage-tracing experiments and compare with Emx1^Cre;Rosa26^YFP embryos. We favoured short-term lineage tracing using Dbx1^LacZ/+ embryos to permanent tracing with Dbx1^Cre, as a delay in recombination was previously shown to prevent labelling of the earliest Dbx1-derived cells (Bielle et al., 2005; Puelles et al., 2016b). We found that a higher proportion of Foxp2⁺ cells derive from the Emx1 lineage (37±5%) than the Dbx1 lineage (14±7%). Conversely, Nurr1⁺ cells are more often Dbx1-lineage derived (43±14%) than Emx1-lineage derived (18±9%). We concluded that Nr4a2⁺ neurons are mostly VP derived, whereas Foxp2⁺ are mostly LP derived, although care should be taken in the interpretation of these data, as (1) β-galactosidase expression is not permanent in Dbx1^LacZ/+ animals and (2) some extent of recombination in the VP is evident in the Emx1^Cre line (sometimes even in subpallium, see Fig. 4D,E) so we cannot formally exclude postmitotic recombination in VP derivatives. Nevertheless, Nr4a2⁺ neurons clearly appear as the most ventral pallial glutamatergic neuron type in our dataset. Our data therefore indicate that neuronal identities in the LP and VP are not sharply distinguished by tracing from Emx1- or Dbx1-expressing progenitor.

Neuronal diversity was proposed to emerge from differences in APs (Telley et al., 2019). To explore AP diversity, we performed a SPRING dimensionality reduction after excluding the cell-cycle associated genes, and subsequently projected cells on the principal curve capturing the main source of transcriptional variability (Fig. 5A). We observed that this axis correlates with the DV position, as illustrated by the sequential distribution of subpallial (Gsx2⁺), VP (Sfrp2⁺) and DP (Emx1⁺) progenitors (Fig. 5B). We found 537 genes showing significant variation along this pseudo-DV axis (Table S6) and used them to discretise APs in seven groups, as we found it was the best trade-off that allowed the correct positioning of the PSB without overclusterising the data (Fig. 5C,D). On the basis of gene expression, pallial APs were found to split into one VP cluster, two LP clusters and one DP cluster. Our approach was validated by the observation that Gsx2, Dlx1/2 or Olig2 were restricted to subpallial clusters, Dbx1 expression was detected almost exclusively in the VP domain (even though its expression is weaker in APs than BPs), and Emx1 was found absent from subpallial clusters and steadily increased in line with the pseudo DV score (Fig. 5D,E, Fig. S5A). Consistently, the expression dynamics of Ascl1, Pax6, Sp8 and Emx2 was found graded along the pseudo DV axis, matching the known expression patterns of these genes (Fig. S5A,B). As expected, the PSB constituted the most remarkable border among APs. We found that it is best defined by the expression of Gsx2 on the subpallial side and Tfap2c or Dmrta2 on the pallial side (Fig. 5E,F, Fig. S5A) consistent with previous findings (Desmaris et al., 2018; Konno et al., 2019). Moreover, we identified multiple additional genes with a sharp PSB expression border, including previously unreported genes such as Gm29260 or Lypd6 (Fig. 5E,F, Fig. S5A,B). However, the overwhelming majority of genes differentially expressed along the pseudo-DV axis did not display sharp borders. For example, Celf4, Sema5a or Fat4 peak at the VP but are also found in adjacent domains (Fig. S5A). Sfrp2, which is often used as a PSB marker, is detected at low levels in all domains, with the few cells apposed to the PSB showing highest expression (Fig. S5A).

We also evaluated whether other sources of variations could account for the diversity of pallial progenitors, especially along the RC axis. We performed PCA and excluded the components correlating with the cell cycle or pseudo-DV axis to project cells on a principal curve (Fig. S5C) that we found reflects the RC axis. This can be exemplified with the patterning genes Pax6, Lhx2 and Nr2f1 (COUP-TFI) for which the inferred DV/RC gradients (Fig. S5C) match those observed in vivo (O'Leary et al., 2007). We identified 159 genes variable along the RC axis, 124 of which were also changing along the DV axis, as illustrated by the expression pattern of Mpped2, Celf4, Fgf15 or Lypd6 (Fig. S5C,D, Table S7). Only a few genes, such as Etv5, were variable strictly in the RC dimension (Fig. S5C), indicating that most of the transcriptional variability among APs is associated with their DV position. Our data therefore provide an atlas of gene expression along the DV axis of the early developing telencephalic VZ and supports the hypothesis that, besides the few genes that sharply delineate the PSB, progenitor identity changes gradually according to the position of cells (Fig. 5G).

The cerebral cortex displays a well-known neurogenic gradient in which neurogenesis is initiated in ventral regions and progressively extends to dorsal territories (Bayer and Altman, 1991). In addition, it was recently shown that neuronal diversity in the developing neocortex arises mostly from the temporal progression of APs (Okamoto et al., 2016; Telley et al., 2019). We therefore queried to what extent the distinction between LP, VP and DP progenitors could mirror that of young versus old dorsal cortex APs.

To test such a possibility, we used the time-series dataset from Telley et al. (2019), which describes the transcriptomic signature of E12 to E15 APs from the DP, and trained a regression model to predict a maturation score for each AP present in our dataset (Fig. 6A). We found a clear difference in the maturation estimate between clusters as the VP displayed the highest predicted score (Fig. 6B). In other words, E12.5 VP progenitors display a more mature transcriptomic profile than their DP counterparts. This also indicates that part of the DV axis transcriptional differences can be explained by the expression of genes responsible for the temporal maturation of dorsal cortex APs.

We found that among the 393 pallial variable genes in our dataset, 102 also vary during dorsal cortex AP temporal maturation, therefore representing a ‘temporal module’ of pallial AP identity. Conversely, the remaining 291 genes form a ‘spatial module’ expected not to show such variation in time (Fig. 6C, Tables S8 and S9). We used genes of the temporal module to investigate whether the temporal-associated signature changes continuously along the spatial and temporal axes. To this aim, we performed correlation between the smoothed gene expression profiles computed in the two datasets and found two blocks of strong correlation between VP/LP.1 and E15 APs on the one hand, and LP.2/DP and E12 APs on the other hand (Fig. 6D). The distinction between LP.1 and LP.2 appears very sharp when considering genes from the temporal module, i.e. genes whose expression is turned on or off when dorsal cortex APs switch from the production of deeper layer to upper layer neurons. We found that within the temporal module, 66 genes are enriched in the E12.5 VP and upregulated at E15 in the DP, as exemplified by the transcription factor Zbtb20 (Fig. 6E,F, Fig. S6D). They can be pictured as a wave of genes expanding dorsally with time. Twenty-eight genes, such as Lrrn1, displayed the opposite dynamics, highly expressed in the DP at E12 and regressing at E15 (Fig. 6E,F, Fig. S6D). The remaining nine genes displayed more complex patterns (Fig. 6E, Table S8).

When focusing on the spatial module (i.e. genes whose expression is constant from E12 to E15), we found far more genes upregulated in the VP than in the DP (225 versus 41, respectively), indicating that the spatial signature of the VP is more pronounced than that of any other pallial domain (Fig. 6G,H, Table S9). Celf4 and Mpped2 exemplify such time-invariant genes enriched either in the VP or the DP, respectively (Fig. 6H). By contrast, only 25 genes were found enriched in the LP.1/LP.2, none of which displayed strong specificity or a sharp pattern of expression (Fig. 6G), indicating that APs in these domains are best defined by the combined low expression of VP- or DP-enriched genes. To further investigate spatio-temporal relationships over different embryonic stages, we took advantage of a recently published dataset (La Manno et al., 2020 preprint) from which we could extract forebrain progenitors collected between E11 and E15.5 (Fig. 6I, Fig. S6A). We found that pallial APs segregate in two groups corresponding to early (E11-E13) and late (E13.5-E15.5) stages, reflecting major changes due to temporal maturation. We computed signature scores corresponding to VP-expressed genes belonging to the spatial or temporal modules, and confirmed that the spatial signature is indeed stable throughout corticogenesis, whereas maturation-associated genes expand in the DP with time (Fig. 6J, Fig. S6C,D). We validated these inferred patterns by performing in situ hybridisation for genes of the spatial and temporal modules at E12.5 and E15.5 (Fig. S6C,D). Altogether, these results indicate that differences between progenitors along the DV axis are explained by the superimposition of a maturation-associated module composed of genes turned on or off between E12 and E15, together with a spatial stage-independent VP or DP identity.

We then attempted to reconstruct developmental trajectories from AP to LN, assuming that most of the glutamatergic neurons we characterised were collected together with their progenitors of origin. Because CRs are known to migrate tangentially over long distances, we excluded them from our lineage analyses.

In order to link neuronal populations to their progenitor pools of origin and identify their lineage-specific intermediate states, we used FateID (Herman et al., 2018), an algorithm that starts from the most differentiated clusters (as defined in Fig. 2D) and progresses backward iteratively to compute a fate probability bias for each cell. From the seven mature populations that were initially seeded, only the two most distant trajectories, leading to Bhlhe22⁺ and Nr4a2⁺ populations, could be confidently reconstructed back to the AP state (Fig. 7A,B, Fig. S7A). For the remaining neuronal types, lineage-specific signatures emerged only at postmitotic stages (Fig. S7B), suggesting that fate acquisition does not occur at the same differentiation step in all lineages or is not strongly transcriptomically encoded. Alternatively, we cannot formally exclude the possibility that some lineage intermediates were missed in the dissection due to cell migration (especially for the Foxp2/Lhx5⁺ population for which the FateID pipeline was unable to reconstruct even the beginning of a trajectory).

Nevertheless, classification scores for APs obtained using this approach supported histological observations, by giving the VP a higher Nr4a2 fate bias while making the DP more likely to generate Bhlhe22⁺ neurons (Fig. 7A,B). Interestingly, we observed a clear distinction between the two lineages within cells at the BP progenitor state (Fig. 7B), with most of those undergoing direct neurogenesis being attributed to the VP/Nr4a2 lineage. To validate this observation, we performed immunostaining for Ki67, Tbr2 and Tbr1 to quantify the amount of cycling cells (Ki67⁺) among BPs (Tbr2⁺/Tbr1⁻). We found the proportion of Ki67⁺ BPs to increase from ventral to dorsal pallial domains (Fig. 7C, Fig. S7C), confirming that lineage amplification through BP cycling is not the same for all cell types.

We identified 231 genes significantly upregulated along the DP trajectory and 426 along the VP trajectory (Fig. 7D, Tables S10 and S11). Those were found at all steps of differentiation, even though signature specificity increased with pseudotime (selected examples are shown in Fig. 7E, Fig. S7D), consistent with the idea that neuronal identity is acquired progressively throughout differentiation. We identified several genes whose expression peaked at the BP stage in the VP trajectory, such as Neurod4, Dbx1, Dleu7 or Svil, but failed to identify counterparts in the DP lineage (Fig. 7E, Fig. S7D). Consistently, in situ hybridisation confirmed that Neurod4 expression is highly enriched in the VP subventricular zone with respect to the DP. In the DP lineage, genes such as Nfib, Neurod2/6 or Pcp4 are sequentially upregulated (Fig. 7E). By contrast, pan-neurogenic genes such as Neurog2, Tubb3 or Mapt (Tau) display identical patterns of expression in the two lineages (Fig. 7E), suggesting that neuronal differentiation proceeds at a similar pace in pseudotime, if not in real time. We then attempted to determine to what extent neuronal fate specification relies on information transmitted from APs or from the onset of lineage-specific gene expression. To do so, we performed a receiver operating characteristic (ROC) curve analysis to identify genes whose expression can best classify BPs between VP and DP lineages (Fig. 7F, Table S12). Out of those best predicting the DP lineage bias, only a few were already differentially expressed among APs such as Emx1. By contrast, the genes best predicting the VP lineage bias in BPs were often found differentially expressed in APs, belonging to both the temporal (e.g. Zbtb20 and Pantr1) and spatial modules (e.g. Epha3 and Meg3) described earlier. These results suggest that commitment to the VP fate mostly relies on genes expressed in APs, with the contribution of both spatial and temporal modules, whereas commitment to the DP fate mostly involves the upregulation of genes that were not already distinguishing APs.

Taken together, our data indicate that significant differences exist between the specification programme of distinct lineages and support the idea that neuronal differentiation proceeds independently of fate acquisition. They also suggest that, at least for some lineages, cortical progenitors are already strongly biased at the AP stage. Our work therefore allows, for the first time, the disentanglement of the sets of genes involved in neuronal differentiation from those controlling neuronal subtypes specification in different regions of the developing cerebral cortex.

In recent years, studies focusing on the establishment of neuronal diversity during cortical development exclusively focused on the neocortex (Loo et al., 2019; Mayer et al., 2018; Telley et al., 2016, 2019; Vitali et al., 2018). Conversely, lateral and ventral regions have been neglected, despite being involved in key physiological processes. The precise contribution of LP and VP progenitors have long been discussed (Medina et al., 2004; Puelles, 2014; Saulnier et al., 2013) and the expression of Nr4a2 has been instrumental in the recent redefinition of LP and VP derivatives (Puelles, 2014; Puelles et al., 2016a). We have shown that most of the early VP-derived neurons (as identified by Dbx1 genetic tracing) correspond to Nr4a2⁺ cells. This comes as a surprise, given that Nr4a2 is only reported to be expressed in the presumptive LP derivatives that are the claustrum and dorsal endopiriform nucleus (Puelles et al., 2016a; Watson and Puelles, 2017), to which the Dbx1 lineage is not supposed to contribute (Puelles et al., 2016b). By contrast, the Fezf1⁺ population we characterised is most likely derived from the LP according to its position in an Emx1-derived domain dorsal to Nr4a2 cells; it also possibly corresponds to neurons that will eventually settle in the lateral amygdala, as expression of Fezf1⁺ has been reported in this region (Hirata et al., 2006; Kurrasch et al., 2007). This population would represent an early-born LP contribution to the amygdala, as previously postulated (Cocas et al., 2009; Gorski et al., 2002). Our findings therefore point to the idea that LP- and VP-derived structures are perhaps not as sharply distinguished as initially anticipated (Puelles et al., 2000).

We captured distinct CR subtypes in our dataset. Because CRs originate from at least four distant and vastly different micro-environments – the septum, hem, VP and thalamic eminence (Bielle et al., 2005; Ruiz-Reig et al., 2017; Takiguchi-Hayashi, 2004) – the concept of CR diversity has long been standing. Although a hallmark of CRs is their absence of Foxg1 expression, we were unable to identify a single gene specifically expressed by all CR subtypes, as exemplified by Reln and Calb2, which are found in the dorsal cortex preplate as well as in some Foxp2 populations. And even though previous work revealed subset-specific features (Griveau et al., 2010; Hanashima et al., 2007), the extent of similarity/differences between CR populations remained largely uncharted. We report that the first and main distinction between CR subtypes corresponds to their lateral or medial origin. Interestingly, VP-derived CRs mostly differ from hem, septum and thalamic eminence CRs by lacking a complete set of genes classically used as CR markers (e.g. Trp73, Lhx1 and Ebf3). By contrast, we could not find genes very specifically expressed by VP-derived CRs, suggesting they represent a default fate, consistent with previous studies indicating that only VP-derived CRs are produced during conditional inactivation of Foxg1 (Hanashima et al., 2007).

Sharp segregation of progenitor domains in the developing spinal cord has long suggested that spatial patterning is a key driver of neuronal diversity. However, in the neocortex, the graded expression of patterning genes, the sequential production of deep and upper layers, and the changes in the transcriptomic signature associated with AP maturation (Okamoto et al., 2016; Telley et al., 2019) clearly point to the importance of temporal mechanisms. If the spatial component has been overlooked in recent studies implementing scRNAseq approaches, one should keep in mind that spatial and temporal regulations are intimately linked. This is perhaps best exemplified by the neurogenic gradient that leads to the generation of ventral neurons earlier than dorsal ones and is associated with opposite gene waves progressing along the DV axis of the brain. The ‘temporal module’ we describe consists of genes whose expression in neocortical progenitors not only changes with time but is also spatially restricted at E12.5 (mostly to the VP). Similarly, the membrane potential of cortical APs, which has been shown to change with temporal maturation, is also different between dorsal and ventral cells at the same developmental stage (Vitali et al., 2018). In the future, it would be interesting to investigate within the DP how much of the variability among APs at a given stage can be attributed to the position cells occupy on the DV axis, and to what extent such differences rely on the temporal maturation gradient. In addition, it remains to be established whether spatial- and temporal-associated gene regulatory networks are independent or cross-regulated. In an evolutionary perspective, a switch from spatial to temporal patterning, involving the co-option of an early VP-specific gene module by DP progenitors, has been proposed to be an important innovation for the acquisition of upper layers in mammals (Luzzati, 2015). Consistently, the loss of function of Zbtb20 or Tfap2c, the expression of which is first restricted to the VP before progressively expanding in dorsal progenitors, result in upper layer neuron generation defects (Pinto et al., 2009; Tonchev et al., 2016). It will thus be interesting in future studies to compare gene expression in cortical progenitors across different species and investigate the extent of conservation, in the composition and dynamics, of the spatial and temporal gene modules.

Along the DV axis of the pallium, APs appear smoothly distributed, with no obvious sharp transition between them, raising the issue of how we should define progenitor domains. It is interesting that a similar question was raised by Wullimann (2017), based on purely neuroanatomical grounds. Perhaps what we currently define as domains do not exist per se but rather APs form a continuum along which the bias towards alternative fates changes according to the DV position and time.

We found that APs occupying the most extreme positions in the gradient (VP and DP) display the strongest spatial signature and bias towards specific terminal fates. By contrast, APs localised at intermediate positions along the DV axis do not show such a strong transcriptional identity that could be transmitted to BPs, explaining why our lineage reconstructions were effective only when comparing the most distant terminal fates in our dataset. However, the near-perfect tangential segregation of Ppp1r14c⁺ and Fezf1⁺ neurons argues in favour of an equivalent spatial segregation of their respective progenitors in the LP. Therefore, we propose that at the cell population level, bias towards a specific lineage progressively decreases, as progenitors are located closer to a transition zone between two domains while the opposite probability towards generating the neighbouring cell type would increase (Fig. 8A). Such a model is supported by the observation that genetic tracing using sharp domain-restricted enhancers results in fuzzy fate maps with graded contribution to adjacent cortical structures (Pattabiraman et al., 2014). At the individual cell level, some progenitors would adopt a deterministic commitment towards one lineage, whereas others, located at transition zones, would undergo a succession of probabilistic choices between two adjacent fates. Clonal fate-tracing techniques provide a powerful experimental paradigm to investigate the behaviour of individual progenitors (Gao et al., 2014; He et al., 2012; Llorca et al., 2019) and we expect future studies taking advantage of such approaches to clarify the link between APs progeny and their position on the DV axis.

The contrast between AP continuity and the discrete nature of the LN also raises the question of the time point at which cell fate is assigned. In the neocortex, graft experiments of DP progenitors indicated that APs display greater plasticity than BPs (Oberst et al., 2019). In our dataset, the bias for Bhlhe22⁺ and Nr4a2⁺ lineages is fairly constant from AP to LN, supporting an early commitment from the AP state. By contrast, the bias for other lineages only emerges at postmitotic stages, suggesting that fate restriction occurs at distinct steps in each lineage, although we cannot formally exclude that commitment indeed occurs at the AP level but remains hidden at the transcriptomic level, as shown in other developmental systems (Weinreb et al., 2020; Ziffra et al., 2020 preprint).

In an attempt to graphically represent our findings, we propose to adapt the classical Waddington metaphor to integrate the influence of both space and time on the same landscape (Fig. 8B). In this representation, AP identity is encoded in its coordinates on a spatio-temporal map. Upon exit from a progenitor state, this coordinate sets the initial position of the cell onto the developing landscape that results from the superimposition of spatial and temporal components. In this model, cells starting from progenitor fields at extremities of the landscape will systematically fall into a basin of attraction that will canalise their differentiation towards the same attractor state (i.e. neuronal fate). By contrast, APs close to the tipping point separating two basins of attraction will be more sensitive to small random fluctuations in their initial state, which might drive them towards either one of two alternative fates. Changes in the regulatory dynamic of genes that underlie the topology of spatial and temporal landscapes result in the modification of the amount and identity of neuronal types being produced from a specific domain. Accordingly, one can also speculate that evolutionary modification of these landscapes might give access to previously unexplored attractors and therefore allow the generation of novel neuronal types.

The following mouse lines were used and maintained on a C57BL/6J background: Dbx1^LacZ (Pierani et al., 2001), Emx1^Cre (Gorski et al., 2002), ΔNp73^Cre (Tissir et al., 2009), Wnt3a^Cre (Yoshida et al., 2006), Rosa26^YFP (Srinivas et al., 2001) and Rosa26^tdTomato (Madisen et al., 2010). All animals were handled in strict accordance with good animal practice, as defined by the national animal welfare bodies, and all mouse work was approved by the French Ministry of Higher Education, Research and Innovation as well as the Animal Experimentation Ethical Committee of Paris Descartes University (CEEA-34, licence number: 18011-2018012612027541).

Four wild-type E12.5 embryos obtained from two distinct litters were collected in ice-cold HBSS. A region spanning from the subpallium to the DP all along the rostro-caudal (RC) axis (see Fig. 1A) was dissected on both hemispheres and dissociated using the Neural Tissue Dissociation Kit (P) (Miltenyi Biotec) and a gentleMACS Octo Dissociator with Heaters following the manufacturer's instructions. Cell clumps and debris were removed by two consecutive rounds of centrifugation for 3 min at 200 g and filtration through 30 µm cell strainers (Miltenyi Biotec). Additional mechanical dissociation and removal of cell aggregates was achieved by gentle pipetting up and down five times using a P200 pipette and Gel Saver tips (QSP), followed by filtration through a 10 µm cell strainer (pluriSelect). Assessment of cell viability and counting were achieved using a MACSQuant flow cytometer (Miltenyi Biotec) in the presence of 7-AAD. The absence of dead cells, debris or doublets in the cell suspension was crosschecked using a haemocytometer and Trypan Blue staining. Approximatively 10,000 cells were loaded on a 10X Genomics Chromium Controller. The entire procedure, from the sacrifice of the pregnant female until loading the controller, was achieved in 2 h. A single-cell barcoded cDNA library was prepared using a Chromium Single Cell 3′ Library and Gel Bead Kit v2, and sequenced at a total depth of 190 million reads on an Illumina NextSeq500 sequencer. Raw sequencing reads were processed to counts matrix with Cell Ranger (version 2.2.1) with default parameters, using the mm10 mouse genome as reference. All statistical analyses were performed under R (version 3.6.3).

To retain only high-quality cells, we applied filters on unique molecular identifier (UMI), gene and mitochondrial counts per cell. We first filtered cells based on the percentage of UMI associated with mitochondrial transcripts and number of detected genes, and retained cells within a three median absolute deviation (MAD) around the population median for these metrics. We further excluded cells having a UMI number above 3 MADs of the population median. Potential doublets were removed using Scrublet (Wolock et al., 2019). Finally, we excluded genes detected in fewer than three cells. Dimensionality reduction was performed using the SPRING tool, a 2D projection using force-directed layout of the k-nearest-neighbour (Weinreb et al., 2018). After dimensionality reduction, 11 cells sharing neighbours between excitatory and inhibitory neurons, and therefore likely corresponding to undetected doublets, were removed manually. The raw counts value for the 4225 high quality cells was library size normalised and scaled using Seurat (v 2.3.4) (Butler et al., 2018).

We reasoned that clustering procedure gives meaningful results only if applied to cell populations in the same broad transcriptional state (e.g. apical progenitors or late neurons). To select cells in the same transcriptional state, we computed apical progenitors (AP), basal progenitors (BP), early neurons (EN) and late neurons (LN) signature scores using lists of published and manually selected genes (Telley et al., 2016; a detailed list is available in the annotated codes, see Data availability section) as input to the Seurat function ‘addModuleScore’. The scores were differentially combined to define and extract the populations of interest by K-means clustering prior to downstream analysis.

To explore neuronal transcriptional diversity, we first applied iterative graph-based clustering using scrattch.hicat (v0.0.16) (Tasic et al., 2016) on glutamatergic and GABAergic neurons separately after removing cell cycle-, ribosomal- or mitochondrial-associated genes from the input count matrix. At each clustering iteration step, we further excluded principal components (PCs) highly correlating with sequencing depth, and with percentage of ribosomal or mitochondrial transcripts detected per cells.

Because variation among some glutamatergic neuron clusters were obviously associated with different maturation states, we further excluded those with a LN score <1. Hierarchical grouping of the remaining clusters was performed with the scrattch.hicat function ‘build_dend’ using the mean gene expression profile across clusters as an input. To define cluster-specific transcriptional signatures, we further subset the clusters core cells as assigned more than 95/100 trials to the same cluster using the centroid classifier of the ‘map_sampling’ function. Top cluster marker genes were identified using a receiver operating characteristic (ROC) analysis implemented in the Seurat function ‘FindMarkers’. Barplot representations of gene expression were built using the scrattch.vis package (v0.0.210).

We first aimed at ordering APs along the main axis of variation by performing a new SPRING dimensionality reduction while excluding cell cycle associated genes. Cells were projected on the principal curve using princurve (v2.1.4) as described previously (Mayer et al., 2018). Genes differentially expressed along this inferred DV axis were identified using the ‘differentialGeneTest’ function of the Monocle package (v2.14) (Trapnell et al., 2014). Progenitor cells were grouped in seven clusters using partition around medoids on the smoothed expression profile of the significantly differentially expressed genes. Genes with similar expression dynamics were grouped by hierarchical clustering.

To investigate sources of variation associated with the RC axis among pallial progenitors, we performed a principal component analysis (PCA) dimensionality reduction using the Seurat function ‘RunPCA’. After excluding the first three PCs correlated (absolute Pearson's cor. >0.25) with the cell cycle or the inferred DV axis, we projected cells along a principal curve spanning the remaining seven PCs and identified significantly variable genes along this new axis using the same approach as described above for the DV axis. Among these 159 genes, we identified those with similar graded expression by fitting a generalised additive model using the ‘gam’ function from the mgcv package (v1.8-33) over both RC and DV axes, and performed PCA on their smoothed expression landscape.

To predict temporal maturation scores for apical progenitors, we took advantage of the Telley et al. (2019) dataset. After retrieving the raw count matrix from the GEO repository (GSE118953), we extracted cells sorted 1 h after FlashTag injection. Based on their expression profile, we excluded medial pallium APs as well as GABAergic interneurons to retain 725 high quality cells (E12, 178 cells; E13, 202 cells; E14, 124 cells; E15, 221 cells). We ran a PCA over all differentially expressed genes between the distinct sampling groups identified with the Wilcoxon rank sum test implemented in the ‘FindAllMarkers’ function. Cells were projected onto a pseudo-maturation axis by fitting a principal curve over the first six PCs. We trained a random-forest regression model implemented in the Caret package (v6.0-84) to predict the pseudo-maturation score (position along the axis) of cells using the scaled expression value of the top 1129 variable genes selected with the ‘FindVariableGenes’ function. The 100 most important features were set as an input to re-optimise the regression model that was used to predict the temporal maturation for APs in our dataset.

To investigate the spatial and temporal components, we first identified genes with variable expression in pallial apical progenitors along the DV axis using the ‘differentialGeneTest’ function as described above, and those differentially expressed between E12 and E15 1 h after FlashTag injection in the Telley et al. (2019) dataset. The ‘temporal module’ was defined as the intersection between the two sets of differentially expressed genes. We then computed the Pearson's correlation coefficient between the smoothed expression matrices over 200 points of temporal genes in our dataset and that of Telley et al. (2019). Genes of the temporal module showing coordinated changes in expression in the two datasets were grouped using hierarchical clustering. The ‘spatial module’ was defined as genes differentially expressed along the DV axis at E12 but that do not show variation between E12 and E15 in the Telley et al. (2019) dataset.

Embryonic brain atlas (La Manno et al., 2020 preprint) expression and metadata were retrieved from http://mousebrain.org/. We used LoomR 0.2.1.9000 (https://github.com/mojaveazure/loomR) to extract cells annotated as ‘Radial glia’ from the nine sampling points ranging from E11 to E15.5. We further retained only those annotated as coming from ‘Forebrain’ tissue. The raw count values were normalised to the library size and scaled using Seurat (v 2.3.4) and UMAP dimensionality reduction was performed using the Seurat function ‘RunUMAP’ on the 30 first PCs. Spatial VP and Temporal VP signatures were computed with the Seurat function ‘addModuleScore’, using genes expressed at the VP and belonging to the ‘spatial’ and ‘temporal’ modules described above.

In order to predict fate bias along the pallial differentiation path, we performed FateID (v0.1.9) (Herman et al., 2018) analysis setting mature pallial neuron clusters (at the exception of CRs) as the ‘attractor’ states. We focused on the two trajectories that were traced back to APs for pseudotime ordering. We restricted our analysis to the confidently assigned cells by selecting those biased towards one lineage with >50% of votes and an absolute difference of votes between the two lineages >25% (Nr4a2 lineage, 617 cells; Bhlhe22 lineage, 520 cells). Cells were projected along a pseudotime axis by fitting a principal curve on the SPRING dimensionality reduction that captures the variation associated with pan-neuronal differentiation. We tested for differential gene expression along pseudotime between the two trajectories using the likelihood ratio test implemented in Monocle's function ‘differentialGeneTest’. Trajectory-specific enrichment was determined by adapting the Area Between Curves (ABC) method from Monocle. Finally, genes with similar smoothed trends were clustered using partition around medoids. To identify genes best predicting early transcriptional differences between the two trajectories at the BP state, we perform ROC analysis as implemented in Seurat's function ‘FindMarkers’ and use the Area Under Curve (AUC) value to rank genes according to their ability to classify BP cells in their respective lineages.

For staging of the embryos, midday of the vaginal plug was considered as embryonic day 0.5 (E0.5). Embryos were collected in ice-cold PBS, dissected and immediately fixed by immersion in 4% paraformaldehyde, 0.12 M phosphate buffer (PB) (pH 7.4) for 2 h at 4°C. Samples were cryoprotected by incubation in 10% sucrose, PB overnight at 4°C, embedded in 7.5% gelatine, 10% sucrose and PB, and frozen by immersion in isopentane cooled at −55°C. 20 µm coronal sections were obtained with a Leica CM3050 cryostat and collected on Superfrost Plus slides (Menzell-Glasser).

The following primary antibodies were used: chick anti-GFP (Aves Labs GFP-1020, 1:1000), chick anti-βGal (Abcam ab9361, 1:1000), rabbit anti-Nfib (Atlas Antibodies HPA003956, 1:1000), rabbit anti-Foxp2 (Abcam ab16046, 1:1000), rabbit anti-Tbr2 (Abcam ab216870, 1:300), rabbit anti-Tbr1 (Abcam ab31940, 1:500), rabbit anti-Ki67 (Abcam ab15580, 1:500), goat anti-Nurr1 (R&D systems AF2156, 1:200), mouse anti-Meis2 (Sigma WH0004212M1, 1:1000) and guinea pig anti-Bhlhb5 (Skaggs et al., 2011, 1:100). The following secondary antibodies were obtained from Jackson ImmunoResearch: donkey anti-chick Alexa-488 (1:1000), donkey anti-mouse Cy5 (1:500), donkey anti-guinea pig Cy5 (1:500), donkey anti-rabbit Cy3 (1:700), donkey anti-rabbit Cy5 (1:500), donkey anti-goat Cy3 (1:700) and donkey anti-goat Cy5 (1:500). Zenon Rabbit IgG labelling kit (Invitrogen) was used to combine several rabbit antibodies according to the manufacturer's instructions. DAPI (1 µg/ml) was used for nuclear staining. Slides were mounted in Vectashield (Vector).

For each gene of interest, a DNA fragment (typically 500-800 bp) was amplified from an embryonic brain cDNA library using Phusion polymerase (Thermo) and a pair of specific primers (Table S1). The promoter sequence of the T7 RNA polymerase (GGTAATACGACTCACTATAGGG) was added in 5′ of the reverse primer. Alternatively, for Dbx1, Emx1, Fgf15, Gsx2, Reln, Rorb, Sfrp2, Sp8 and Trp73, a plasmid containing part of the cDNA was linearised by enzymatic restriction. Antisense digoxigenin-labelled RNA probes were then obtained by in vitro transcription using T7, T3 or SP6 RNA polymerase (New England Biolabs) and digRNA labelling mix (Roche). In situ hybridisation was carried out as previously described (Schaeren-Wiemers and Gerfin-Moser, 1993) using a hybridisation buffer composed of 50% formamide, 5×SSC, 1×Denhardt's, 10% dextran sulfate, 0.5 mg/ml yeast RNA and 0.25 mg/ml herring sperm DNA. Probes were detected using an anti-digoxigenin antibody coupled to alkaline phosphatase (Roche) and NBT/BCIP (Roche) as substrates. Slides were mounted in Mowiol.

In situ hybridisation images were obtained using a Hamamatsu Nanozoomer 2.0 slide scanner with a 20× objective. Immunofluorescence images were acquired using a Leica SP8 confocal microscope with a 40× objective. Quantifications were achieved by counting the proportion of GFP⁺ cells among the marker-positive population using the ImageJ software. Data are mean±s.d., individual values are plotted.

The authors acknowledge Aurelia Dujardin, Camelia Piat and the Animalliance platform of the Imagine Institute for animal care; the SFR Necker imaging and histology platforms for help with acquisition; Fanny Coulpier from the IBENS genomics platform for sequencing; Christine Bole and Aurore Pouliet from the Imagine Institute genomics platform for helpful discussions; Frédéric Tores, Jean-Marc Plaza and Patrick Nitschke from the Imagine Institute bioinformatics platform for implementing the offline version of SPRING Viewer; 10X Genomics and Miltenyi Biotec technical support for their help in the implementation of scRNAseq. We thank all members of the Pierani lab, Marine Luka and members of the Ménager lab, as well as members of the Rausell lab for helpful discussions. We are grateful to Anne Teissier and Oriane Blanquie for critical reading of the manuscript.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.197962