Radial glia (RG) in the neocortex sequentially generate distinct subtypes of projection neurons, accounting for the diversity and complex assembly of cortical neural circuits. Mechanisms that drive the rapid and precise temporal progression of RG are beginning to be elucidated. Here, we reveal that the RG-specific transcriptional regulator PRDM16 promotes the transition of early to late phase of neurogenesis in the mouse neocortex. Loss of Prdm16 delays the timely progression of RG, leading to defective cortical laminar organization. Our genomic analyses demonstrate that PRDM16 regulates a subset of genes that are dynamically expressed between early and late neurogenesis. We show that PRDM16 suppresses target gene expression through limiting chromatin accessibility of permissive enhancers. We further confirm that crucial target genes regulated by PRDM16 are neuronal specification genes, cell cycle regulators and molecules required for neuronal migration. These findings provide evidence to support the finding that neural progenitors temporally shift the gene expression program to achieve neural cell diversity.
Radial glia (RG) in the developing mammalian cerebral cortex are neural stem cells that give rise to all excitatory neurons (Anthony et al., 2004; Kriegstein and Alvarez-Buylla, 2009). During the peak phase of neurogenesis, RG divide asymmetrically to produce a self-renewing RG and a neuron or a transit-amplifying intermediate progenitor (IP) that divides again to produce more neurons (Noctor et al., 2004). On each embryonic day (E), starting at E11.5, RG generate a new laminar layer with distinct neuronal subtypes. Layer (L) 6 neurons are born first at E12.5, followed by L5 (E13.5), L4 (E14.5) and L2/3 (E15.5) neurons. The newborn neurons migrate along the radial fiber of the mother RG, passing through and positioning on top of earlier-born neurons (Kwan et al., 2012; Angevine and Sidman, 1961; Okano and Temple, 2009). Thus, the identity and laminar position of neuronal subtypes are tightly linked to their birthdates.
The competence of progenitors becomes progressively restricted (Desai and McConnell, 2000; Gaspard et al., 2008; Gao et al., 2014; Frantz and McConnell, 1996), and the temporal identity of neural progenitors is controlled by both extrinsic and intrinsic mechanisms (McConnell and Kaznowski, 1991; Ge et al., 2006; Dennis et al., 2017; Fukumitsu et al., 2006; Shen et al., 2006; Hsu et al., 2015; Chenn and Walsh, 2002; Yoon et al., 2017; Vitali et al., 2018). Recent single cell transcriptome analyses revealed that RG at different neurogenesis stages may possess temporal gene expression programs (Telley et al., 2019; Yuzwa et al., 2017). However, the mechanisms of how such temporal expression is established remain to be explored.
In differentiating neurons, cell-specific transcription factors and their regulated transcriptional cascades further guide neuronal specification, migration and circuit assembly (Greig et al., 2013). For example, complex interplay between TBR1, FEZF2 and SATB2 guide specification of corticothalamic, subcerebral and callosal neurons in deep-, mid- and upper-cortical layers (McKenna et al., 2015; Srinivasan et al., 2012). Two related POU domain transcription factors, POU3F2 (BRN2) and POU3F3 (BRN1), are required for determining the identity and migration of upper layer neurons (McEvilly et al., 2002; Sugitani et al., 2002). Although key neuronal specification factors are used as subtype-specific markers (Molyneaux et al., 2007), the mRNAs of these factors already appear in neural progenitors (Yoon et al., 2017; Zahr et al., 2018). Whether transcription of neuronal specification genes in the RG is important for rapid daughter cell differentiation after asymmetric division remains to be demonstrated.
The choroid plexus (ChP) protects the brain via the blood-cerebrospinal fluid (CSF) barrier and regulates CSF composition via specific transfer and secretion (Johansson, 2014; Lehtinen et al., 2013). Signaling molecules in CSF (e.g. SHH, IGF1, WNT4, TGM2 and FGF2) are delivered to NSCs and influence NSC behavior (Johansson, 2014; Imayoshi et al., 2008; Lehtinen et al., 2011; Johansson et al., 2013). Mechanisms and factors controlling development of the ChP are not fully understood.
The PR domain-containing (PRDM) family protein PRDM16 is a key transcriptional regulator in diverse cell types (Aguilo et al., 2011; Chuikov et al., 2010; Kajimura et al., 2008). In brain development and homeostasis, PRDM16 was shown to control neural stem cell maintenance and proliferation (Chuikov et al., 2010; Shimada et al., 2017; Chui et al., 2020), IP proliferation (Baizabal et al., 2018), neuronal migration (Inoue et al., 2017; Baizabal et al., 2018; Chui et al., 2020) and ependymal cell differentiation (Shimada et al., 2017). PRDM16 regulates genes involved in reactive oxygen species (ROS) levels (Chuikov et al., 2010; Inoue et al., 2017; Chui et al., 2020) and epigenetic states of its bound enhancers (Baizabal et al., 2018). The PRDM16 protein (Fig. S1A) contains a PR domain that possesses intrinsic Histone H3K4 (Zhou et al., 2016) and H3K9 methyltransferase activity (Pinheiro et al., 2012), two potential DNA binding zinc-finger clusters (Nishikata et al., 2003) and interaction motifs for the co-repressors CtBP1/2. The transcriptional activity of PRDM16 is context-dependent (Chi and Cohen, 2016), as it activates gene expression when associated with activators and represses gene expression when interacting with co-repressors.
In this study, we investigated the mechanisms by which PRDM16 regulates RG property in the developing mouse brain. We show that Prdm16 deficiency compromises the development of ChP and cortical layer lamination and that specific depletion of Prdm16 from the forebrain resulted in similar defects in the neocortex. Although the number of RG in knockout animals appeared to be normal, there was a delay in the progression of RG from producing earlier-born L5 neurons to later-born L2-4 neurons. Our genomic data demonstrated that PRDM16 primarily represses target gene expression through modulating permissive enhancers. By integrating public scRNA-seq data with our RNA-seq and ChIP-seq analysis, we identified a subset of PRDM16-regulated genes that are dynamically expressed from early to late neurogenesis. We propose a model in which PRDM16 helps to set up the temporal transcriptional landscape of RG. As a proof of principle, we confirmed that two crucial target genes of PRDM16, Cdkn1c and Flrt3, the expression of which declines in late neurogenesis, are responsible for the regulation of progenitor proliferation and neuronal migration, respectively.
PRDM16 is required for neocortical development and ChP formation
To assess the function of Prdm16 in the developing brain, we used three multifunctional conditional gene trap (cGT) alleles (Strassman et al., 2017) (Fig. S1B). The Prdm16cGT (cGT) and Prdm16cGTreinv (cGTreinv) mouse strains produce a null allele and will be referred as Prdm16 KO mutants. To examine PRDM16 activity in the neocortex, we depleted Prdm16 expression in the forebrain using the Emx1tm1(cre)Krj/J (Emx1IREScre) deleter strain (Gorski et al., 2002) and the conditional Prdm16cGTinv (cGTinv) strain. cGTinv will be referred to as Prdm16 cKO throughout this manuscript. The Prdm16 transcript is detectable in E9.5 brain (Fig. S1C; (Horn et al., 2011). At E13.5, PRDM16 has specific expression in the ChP and in the ventricular zone (VZ), where it colocalizes with the RG marker SOX2 (Fig. S1D). In KO animals, PRDM16 staining is lost in the entire brain (Fig. S1D), whereas in cKO mutants it is depleted in the dorsal telencephalon but remains expressed in the ventral telencephalon and the ChP (Fig. S1E).
We first analyzed the cortical laminar organization of the KO brains by labeling cortical neurons with SATB2 for the upper-layer (L2-4, II-IV) and CTIP2 (also known as BCL11B) for the mid-layer (L5, V; strong CTIP2) and the deep-layer (L6, VI; weak CTIP2). At postnatal day (P) 0, mutant cortices showed expansion of the CTIP2+ layer, accompanied by thinning of the SATB2+ upper layer (Fig. 1A,B), compared with control cortices. Some SATB2+ neurons scattered inside the deep layer, indicative of failed migration. Similarly, at E15.5, when upper-layer neurons were born, the number of SATB2+ neurons decreased and the mid-layer neurons labeled with FEZF2 and CTIP2 expanded in the mutant (Fig. S1F,G). We confirmed the change using reverse transcription followed by quantitative PCR (RT-qPCR). The levels of the two L5 genes increased to ∼150%, whereas those of the L2-4 genes decreased to 50-70% (Fig. 1C), indicating that gain of mid-layer neurons roughly compensates for loss of upper-layer neurons at E15.5. Hence, the Prdm16 KO cortex displays two types of defects: overproduction of mid-layer neurons, and reduced production and defective migration of upper layer neurons.
In Prdm16 KO brains, the prospective ChP in the lateral and the third ventricles are dramatically reduced (Fig. S1D,H,I; Bjork et al., 2010b; Strassman et al., 2017), pointing to an essential role of PRDM16 in ChP development. Together, the phenotypic analyses in Prdm16 KO mutant suggest that PRDM16 controls brain development in at least two brain areas, the neocortex and the ChP.
Expression of Prdm16 in the forebrain regulates cortical laminar organization
To test a direct role of PRDM16 in cortical development, we analyzed Emx1IREScre-mediated Prdm16 cKO mutants, in which Prdm16 is depleted in the forebrain (Fig. S1E). The Prdm16 cKO animals survive to adulthood, allowing examination of postnatal stages. At E15.5, Prdm16 cKO cortices displayed defects similar to Prdm16 KO embryos: an increase in the number of CTIP2+ and FEZF2+ neurons and a reduction in the number of SATB2+ neurons (Fig. 1D,E). At P15, the cKO cortex showed similar changes (Fig. 1F-H), and the TBR1-labeled deep-layer is unaffected. A subset of SATB2+ neurons failed to migrate and remained below the cortex, as a chunk of grey matter cells (heterotopia) (Fig. 1F,G). Thus, the forebrain depletion of Prdm16 led to the same effects on cortical laminar organization as the null KO did: reciprocal changes of L5 and L2-4 neurons and failure of upper-layer neuron migration. This result confirms that the laminar organization phenotypes in the mutant cortex are due to loss of Prdm16 in the forebrain.
PRDM16 regulates the transition of mid-to-late neurogenesis
We sought to understand the causes of Prdm16 mutant phenotypes. Given that PRDM16 is an RG-specific factor, PRDM16 may control neurogenesis through modulating intrinsic properties of RG. We reasoned that two possibilities could lead to the increase of L5 neurons and decrease of L2-4 neurons. First, if loss of Prdm16 delayed the transition of neurogenesis from E13.5 to E14.5, mutant RG would produce L5 neurons even after E13.5, which could result in fewer L2-4 neurons. Second, if loss of Prdm16 increased proliferation of RG at E13.5 and reduced it at E14.5 and later, more mid-layer neurons could be produced at E13.5 and fewer upper-layer neurons produced at a later time. We attempted to test these possibilities.
To test whether PRDM16 controls the timing of RG transition, we traced RG daughter cell fate by injecting pregnant mice with 5-bromo-2'-deoxyuridine (BrdU) at E14.5 and 5-ethynyl-2′-deoxyuridine (EdU) at E15.5, and examined the distribution and cell identity of BrdU and EdU cells in the P5 cortex. CTIP2+ L5 neurons are born at E13.5 and should not be labeled with BrdU or EdU (Desai and McConnell, 2000; Gaspard et al., 2008). This was seen in the control cortex (Fig. 2A,B). In contrast, supernumerary CTIP2+BrdU+ neurons appeared in the mutant cortex, suggesting CTIP2+ neurons were produced at or after E14.5 in the mutant. Notably, the number of the CTIP2+BrdU− cells did not change in the mutant, indicating that production of CTIP2+ neurons before E14.5 was normal. These results demonstrate that some of Prdm16 mutant RG failed to transit from E13.5 to E14.5 and continued to produce CTIP2+ neurons at E14.5 (Fig. 2C).
BrdU+ and EdU+ cells were found in the heterotopia and the deep layer (Fig. S2A-C), confirming failed migration of some of the upper-layer neurons. None of the CTIP2+BrdU+ cells were retained in the deep layer (Fig. 2A) or the heterotopia (Fig. S2B), suggesting that even the latter-produced CTIP2+ neurons migrate normally and that the migration failure is specific to upper-layer neurons.
PRDM16 promotes indirect neurogenesis and production of IPs
Next, we checked the progenitor populations in Prdm16 mutant. We examined RG and IP cell counts at E15.5 by labeling RG with PAX6 and IPs with TBR2 (EOMES). There was a significant reduction in the number of IPs in the cKO cortex, whereas PAX6+ RG were not affected (Fig. 3A,B), suggesting that PRDM16 is required for IP cell production.
RG divide and produce neurons directly (direct neurogenesis) or indirectly through IP cells (indirect neurogenesis). Although both division modes are present throughout cortical neurogenesis, indirect neurogenesis becomes more predominant in later stages (Yuzwa et al., 2017) (Vitali et al., 2018). As the IP cell number is reduced in the Prdm16 mutant, we examined whether PRDM16 promotes indirect neurogenesis. We electroporated a short-hairpin RNA (shRNA) construct targeting Prdm16 (Bjork et al., 2010a) and an empty shRNA construct to wild-type E15.5 brains in utero. shPrdm16-transfected cells led to significantly fewer TBR2+ cells compared with those transfected with the shRNA vector (Fig. 3C,D). These results suggest that PRDM16 promotes indirect neurogenesis.
The duration of S phase to M phase in E15 RG is 6 h (Takahashi et al., 1995). If a daughter cell of RG has exited the cell cycle, this cell will be labeled with EdU but not with the cell cycle marker KI67 (MKI67) after 6 h. EdU+KI67− cells will reflect the collection of post-mitotic neurons generated from direct neurogenesis and cells that do not re-enter the cell cycle. We injected EdU into pregnant mice with embryos at E15 and analyzed the brains after 12 h. Compared with the control, the Prdm16 cKO cortex displayed a significant increase in the percentage of EdU+KI67− cells among EdU+ cells in the VZ/subventricular zone (SVZ) region, indicating more cells exiting the cell cycle in the mutant (Fig. 3E,F).
Moreover, we detected a significant decrease of KI67+ cells specifically in the mutant SVZ (Fig. 3G), suggesting that PRDM16 also regulates IP cell proliferation. To confirm this, we injected animals with EdU at E15.5 and collected brains 2 h after. The fraction of EdU+TBR2+ cells over TBR2+ cells is significantly less in the mutant compared with the control (Fig. 3H,I; Fig. S3), showing fewer TBR2+ cells in S phase in the mutant. In addition, mitotic cells labeled with phosphorylated Histone 3 (PH3) were reduced in the basal domain of the ventricle in mutant cortices (Fig. 3J,K). This evidence indicate that loss of Prdm16 affects IP cell proliferation.
We next examined cell counts and proliferation of RG and IP cells at E13.5. In contrast to E15.5, neither the PAX6+ nor the TBR2+ cells showed a significant change in cKO cortices (Fig. S4A,B). Staining with PH3 confirmed no change in the number of mitotic cells (Fig. S4C,D). EdU labeling from E13 to E13.5 showed that the fraction of EdU+KI67− cells over all EdU+ cells in cKO cortices appeared to be similar to the control (Fig. S4E,F), indicating that cells exit the cell cycle normally in the mutant at this stage. Another RG marker, SOX2, did not show any change (Fig. S4G). Together, these results suggest that PRDM16 regulates RG neurogenesis in a stage-specific manner: first, it promotes the temporal transition of RG between E13.5 and E14.5; second, it promotes indirect neurogenesis and IP proliferation during late-neurogenesis.
We also attempted to check the effects by overexpressing PRDM16. Surprisingly, overexpressing PRDM16 resulted in severe migration defects and fewer TBR2+ IP cells, compared with the control injection (Fig. S4H-J). This means that the level of PRDM16 is important for regulating neuronal migration and IP cell production.
PRDM16 modulates levels of neuronal specification genes in RG
PRDM16 is a transcription factor specifically expressed in the RG in the neocortex. Because the cell fate change in Prdm16 mutants occurs at the transition between mid-layer neuron production and upper-layer neuron production (Fig. 2), a direct regulatory role of PRDM16 on RG gene expression may be crucial at E13.5. To test this, we generated RNA-seq data from E13.5 control and KO mutant forebrains (FB) (Fig. 4A), and identified 35 downregulated and 47 upregulated genes in KO versus control (cutoff: P-value<0.05 and fold-change>1.5-fold; Fig. S5A). A recent study published the RNA-seq data of sorted RG, IPs and cortical neurons (CNs) from E15.5 control and Emx1IREScre-mediated Prdm16 cKO cortices (Baizabal et al., 2018). We confirmed that most of the deregulated genes in the mutant FB were deregulated in the E15.5 mutant RG (Fig. S5A,B). Consistent with RG-specific expression of Prdm16, more genes changed expression level in mutant RG than those in IPs and neurons.
As the Prdm16 mutant cortex displays abnormal laminar layer organization, we next examined expression of layer marker genes (Molyneaux et al., 2007; Zahr et al., 2018) in the RNA-seq data (Fig. 4A). Several upper-layer genes showed decreased expression in the E13.5 KO FB (Fig. 4B). Using a Limma-based gene set testing (see Materials and Methods), we confirmed a significant downregulation of the upper-layer markers as a gene set (Pdownregulation=0.0007). As the upper-layer neurons are not specified at E13.5, the expression change of upper-layer genes likely occurred in progenitor cells. In contrast, neither the mid-layer nor the deep-layer genes displayed significant change as a group in the KO FB, although the expression of three genes, Pcp4, Otx1 and Fezf2, showed mild increase (Fig. 4B). We also analyzed changes of the layer-specification genes in the published E15.5 RNA-seq data (Baizabal et al., 2018) and observed a similar trend: the upper-layer genes are significantly downregulated (Fig. S5B) and the mid-layer genes as a group showed significant upregulation in mutant RG and IPs, but not in CN, consistent with our hypothesis that these genes are regulated by PRDM16 in progenitor cells.
At E13.5 the mid- and deep-layer neurons are specified, and many of the mid-layer genes are also expressed in the deep-layer. These may mask cell-type specific expression change for the mid- and deep-layer genes from whole FB RNA-seq data. To verify cell-specific changes, we examined expression of Pou3f2 and Unc5d, two upper-layer genes, and Pcp4 and Fezf2, two mid-layer markers in Emx1-Cre Prdm16 cKO mutants using RNA in situ hybridization (Fig. 4C,D). The signal intensity of Pou3f2 and Unc5d was reduced in the mutant VZ/SVZ area at both stages (Fig. S5C), confirming the reduction of their expression in mutant progenitors. Expression of the mid-layer gene Pcp4 showed a mild but significant increase in the mutant VZ/SVZ at E15.5. Fezf2 displayed a trend of increased expression but is shy of significance. This suggests that derepression of Fezf2 may be limited in a small population of cells. These results demonstrate that neuronal genes are already expressed in progenitors and that the normal expression of these genes is disrupted by Prdm16 depletion.
PRDM16 mainly functions as a transcriptional repressor in RG
To determine direct targets of PRDM16 in the developing brain and investigate how the targets are regulated, we performed chromatin immunoprecipitation followed by deep sequencing (ChIP-seq) from E13.5 heads of wild-type and Prdm16 KO animals. Using an Irreproducibility Discovery Rate (IDR) pipeline (see Materials and Methods), we identified 2319 confident peaks (IDR<5%), of which 42.9% were mapped to intergenic regions, 32.7% to introns and only 18.3% close to promoters. This indicates that PRDM16 mainly targets distal enhancers (Fig. 5A). Gene Ontology (GO) analysis shows that PRDM16-bound genes are enriched for nervous system development, migration signaling and RG function (Fig. S6A,B). We also compared our E13.5 whole-head ChIP with the published E15.5 cortex ChIP data (Baizabal et al., 2018), and found that ∼30% (760) of the E13.5 peaks overlap with the E15.5 peaks (Fig. S6C). The overlapping sites represent continuous binding by PRDM16 between E13.5 and E15.5.
Transcriptional activity of PRDM16 is context-dependent (Chi and Cohen, 2016). We asked how the targets in the neocortex are regulated by PRDM16. To address this, we examined chromatin accessibility of PRDM16-targeted enhancers from E13.5 control and mutant cortices by using Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) (Buenrostro et al., 2013). High ATAC-seq signal in the genome correlates with active cis-regulatory elements (Daugherty et al., 2017). To examine accessibility of PRDM16-bound regions, we plotted mean reads per million (RPM) of PRDM16 ChIP peaks in a metaplot and ranked the peaks in a heatmap, and then compared the ATAC-seq signal intensity in the wild type and mutant for each peak (Fig. 5B). The PRDM16-bound regions display high ATAC-seq intensity, suggesting that these regions are regulatory sites. Interestingly, the ATAC-seq signal at these sites became higher in the Prdm16 mutant, suggesting that loss of Prdm16 led to increased chromatin accessibility at PRDM16 binding sites. We quantified changes of ATAC-seq coverage on the PRDM16 ChIP-seq peaks between control and mutant and found that 226 and 189 peaks showed increased and reduced coverage, respectively [false discovery rate (FDR)<0.2 and fold-change>1.4-fold] (Fig. 5C,D). We assessed expression changes for the gene loci that associate with accessibility changes by combining the ATAC-seq data with our E13.5 FB RNA-seq data (Fig. 5C) or the published E15.5 RG RNA-seq data (Fig. 5D). In both cases, the majority of the significantly upregulated genes in the mutant displayed increased chromatin accessibility at PRDM16 binding sites, whereas downregulated ones do not show such a trend.
We then applied a gene set testing for all PRDM16 targets as a set, and found that the targets at both E13.5 and E15.5 have a significant trend of derepression in mutant RG (P<0.001) and IPs (P<0.001), but not in CNs (P=0.19) (Fig. S6D), confirming that many targets are normally repressed by PRDM16. We looked at a group of genes called RG core identity genes (Yuzwa et al., 2017) that are highly expressed in RG throughout neurogenesis. We found that 20 out of the 90 RG identity genes showed deregulation in the Prdm16 mutant FB. Only the subset bound by PRDM16 became upregulated in the mutant FB or the mutant RG and IPs (Fig. S6E,F). Validation on one of the genes, Veph1, by RT-qPCR and in situ hybridization confirmed derepression of Veph1 in Prdm16 mutants (Fig. S6G-I). We therefore conclude that PRDM16 primarily acts as a repressor in RG through limiting chromatin accessibility.
PRDM16 directly represses mid-layer genes including Fezf2
As PRDM16 represses transcription, expression of its targets in RG should be relatively low. We asked in which cell types the target genes have higher expression. To address this, we re-analyzed the published scRNA-seq data from E13.5 (Yuzwa et al., 2017) to obtain cell-type specific transcriptomes. We identified six cell-type-specific clusters and assigned cell type to each cluster based on the presence of known markers (Fig. S7A,B). Consistent with the previous finding (Zahr et al., 2018), the RG and the IP clusters express many neuronal marker genes (Fig. S7B). We then plotted the percentage of cells that contain high expression of PRDM16 targets per cell in each cluster (cutoff: total normalized reads per cell >200, see Materials and Methods; Fig. 6A; Fig. S7C). In addition to the RG cluster, the mid- and deep-layer neuron clusters showed the highest enrichment of PRDM16 targets, suggesting that many PRDM16 targets are highly expressed in mid- and deep-layer neurons.
One of the PRDM16 targets is Fezf2 (Fig. 6B), a mid-layer neuron determinant (Molyneaux et al., 2005). PRDM16 binds to the distal enhancer of Fezf2 known to drive Fezf2 expression in RG (Shim et al., 2012). We confirmed that Fezf2 is de-repressed in Prdm16 KO mutant neural stem cells, using primary NSC culture (Fig. 6C). We then tested responsiveness of this RG enhancer to PRDM16 using a luciferase reporter (Fig. 6D). PRDM16 fused with the VP64 activation domain induced higher expression of the Fezf2 reporter than the pGL3-promoter, confirming PRDM16 binding to this enhancer. However, PRDM16 (PRDM16-FL) or the truncated version (PRDM16-PRdeletion, lack of the PR domain) did not regulate the reporter (Fig. 6D). We reasoned that PRDM16 may require chromatin context for its regulatory activity which is lacking in transient transfection assay. To overcome this, we measured endogenous level of Fezf2 mRNA by RT-qPCR from the N2A cell lines that express the control construct, PRDM16-FL or PRDM16-PRdeletion (Fig. 6E). PRDM16-FL could suppress Fezf2 expression, suggesting that PRDM16 needs endogenous chromatin context or other cis-regulatory element(s) to repress Fezf2. Expression of Fezf2 is not changed in cells expressing PRDM16-PRdeletion, indicating that the PR domain is essential for its repressive activity.
PRDM16 influences temporal dynamics of gene expression in RG
RG produce mid-layer neurons at E13.5 and generate upper-layer neurons from E14.5 to E16.5. We hypothesized that the gene expression program of E13.5 RG differs from that of E15.5 RG to confer temporal identity. As Prdm16 mutants show a delay of RG transition from mid- to late-neurogenesis, PRDM16 may influence the dynamics of gene expression in RG. To test this, we determined dynamically expressed genes in RG between E13.5 and E15.5, by performing differential gene expression analysis between the E13.5 and the E15.5 RG clusters from the published scRNA-seq data (see Materials and Methods). We found that 120 and 248 genes showed higher expression at E13.5 and at E15.5, respectively (FDR<0.2, FC>1.4-fold) (Fig. 7A). We then examined expression changes of these genes in Prdm16 mutant RG [Prdm16 cKO RNA-seq data from Baizabal et al. (2018)], and noticed that 25 of them showed significant changes (FDR<0.05) (Fig. 7B; Fig. S8A). Among these genes, the upregulated genes all contain PRDM16 ChIP peaks, consistent with the repressive role of PRDM16 on gene expression.
In the Prdm16 mutant cortex, the transition to late neurogenesis is disrupted and many of the upper-layer neurons fail to migrate. Thus, we focused on the genes that display higher expression at E13.5 than E15.5. Among three of these genes, Cdkn1c encodes the cell cycle regulator p57KIP2 that suppresses progenitor proliferation in early neurogenesis (Mairet-Coello et al., 2012). Cdkn1c is more strongly expressed in IP than in RG (Fig. S8B,C), suggesting that downregulation of Cdkn1c in later stages may be required to allow higher proliferation of IPs. Another target gene Flrt3 encodes fibronectin leucine rich transmembrane protein 3, a repulsive cue for the UNC5 family receptors in guiding cell migration (Yamagishi et al., 2011). Expression of Flrt3 is also significantly upregulated in Prdm16 mutant RG (Fig. S8D). As Unc5d is expressed in upper-layer neurons, the action between Unc5d and Flrt3 provides a possible mechanism specific for upper-layer neuron migration.
To validate the prediction from the bioinformatic analysis, we examined expression of Cdkn1c and Flrt3 in the control and cKO cortices by RNA in situ hybridization. Both genes displayed higher expression in RG at E13.5 than E15.5 (Fig. 7C; quantified in Fig. S8E). Expression of Cdkn1c is undetectable in the control cortex at E15.5, but became stronger at the VZ/SVZ area in the mutant. Similarly, expression of Flrt3 increased in the mutant VZ/SVZ compared with the control. We then screened the Cdkn1c and the Flrt3 gene loci for candidate enhancers that can be directly bound by PRDM16. Using a luciferase reporter assay in N2A cells, we identified two distal enhancers for Cdkn1c and one enhancer for Flrt3 that responded to ectopically expressed PRDM16-VP64 fusion protein (Fig. 7D-F). Thus, PRDM16 is able to bind regulatory regions of the predicted target genes.
To test whether Cdkn1c and Flrt3 are important for PRDM16 regulation on progenitor proliferation and neuronal migration, respectively, we examined their genetic interactions with Prdm16. We first designed shRNA constructs to specifically knockdown Cdkn1c or Flrt3, and confirmed the knockdown efficiency of these constructs in cultured cells (Fig. S8F-G). We then electroporated shPrdm16 alone and shPrdm16 in combination with shCdkn1c or shFlrt3 into wild-type E14.5 brains in utero. To examine the effects on IP proliferation, we also injected EdU 2 h before harvesting the E16.5 brains. Similar to reduced IP proliferation in Prdm16 cKO mutants, shPrdm16 transfection decreased the fraction of GFP+TBR2+EdU+ cells in GFP+TBR2+ cells compared with the control transfection (Fig. 7G,H; Fig. S9). Additional knockdown of Cdkn1c could partly restore proliferative IPs, suggesting that Cdkn1c is an important target gene responsible for progenitor proliferation.
To score migration defects, we harvested P5 brains and analyzed the distribution of the transfected cells labeled with GFP. As expected, shPrdm16 treatment led to upper-layer neurons remaining in the white matter (Fig. 7I,J; Fig. S8H), whereas shFlrt3 co-injection significantly suppressed this effect. This result confirms that PRDM16 controls upper-layer neuron migration at least partly through regulating expression of Flrt3. Together, our bioinformatic analyses and functional validation demonstrate that PRDM16 regulates expression of a subset of temporally-dynamic genes that mediate its role in promoting temporal transition of RG (Fig. 8).
Our results demonstrate that PRDM16 is a crucial transcriptional regulator that controls the gene expression program of RG during cortical neurogenesis. To ensure the timed progression of RG between early and late phases of neurogenesis, PRDM16 regulates a stage-specific gene expression program that includes at least three classes of target genes (Fig. 8): (1) mid-layer determinants to allow timely upregulation of upper-layer genes in RG; (2) cell cycle inhibitors to allow higher proliferation of IPs at later neurogenesis; (3) genes encoding guidance cues needed for upper-layer neuronal migration.
Baizabal et al. recently reported that PRDM16 regulates upper-layer neuron production and migration, but found no effect on mid-layer neurons (Baizabal et al., 2018). However, we demonstrated that there is a prolonged production of mid-layer neurons in the Prdm16 cKO cortex in addition to its effects on upper-layer neurons. The discrepancy may result from the methods used to assess cell fates. When assessing neuronal population in the cortex, we distinguished L5 mid-layer neurons from L6 deep-layer neurons, whereas Baizabal et al. assessed L5 and L6 neurons as a whole population. Our RT-qPCR and in situ hybridization assays provide additional evidence showing the reciprocal changes of mid-layer and upper-layer marker genes. How RG change their temporal identity is an interesting question. We found that a subset of RG-specific genes show dynamic expression between mid- to late- neurogenesis and demonstrate that timely control of gene expression by PRDM16 is important for temporal progression to late neurogenesis.
Temporal progression of RG regulated by PRDM16
The understanding of temporal control of RG has been augmented over the years. A number of transcription factors and epigenetic regulators have been shown to control the timing of cortical neuronal specification (McConnell and Kaznowski, 1991; Ge et al., 2006; Dennis et al., 2017; Fukumitsu et al., 2006; Shen et al., 2006; Hsu et al., 2015; Chenn and Walsh, 2002). Progressive hyperpolarization of the cell membrane of RG also regulates sequential generation of neuronal subtypes through modulating Wnt signaling (Vitali et al., 2018). Post-transcriptional mechanisms, including translational repression (Zahr et al., 2018) and the N6-methyladenosine (m6A) RNA modification (Yoon et al., 2017), promote RG progression through preventing precocious production of neuronal proteins. However, many questions remain. For example, the gene expression program in RG is primed for daughter cell genes, but how is the priming status of RG established? How does the pre-established transcriptional program impact daughter cell fate? We found that PRDM16 regulates two aspects of RG temporal identity, the balance of direct and indirect neurogenesis and the transition from mid- to late-neurogenesis. Prdm16 mutant RG display abnormal expression of mid- and upper-layer genes and several temporally-dynamic genes involved in proliferation (Id2, Cdkn1c) and migration (Flrt3, Dcx, Sparcl1). These results suggest that PRDM16 is involved in setting up the primed gene expression program of RG.
PRDM16 is expressed throughout cortical neurogenesis. An important question is how its activity is triggered at the onset of mid-to-late neurogenesis transition. Interestingly, Prdm16 co-clusters with Slc1a3, a regulator of metabolism of glutamate and ion flux (Vandenberg and Ryan, 2013). It will be of interest to test a potential functional link between SLC1A3 and PRDM16 in integrating extrinsic and intrinsic signals. Alternatively, PRDM16 uses other mechanisms to sense environmental cues. One such candidate is the level of ROS. It has been reported that PRDM16 integrates ROS signaling with stem cell behavior (Inoue et al., 2017; Chui et al., 2020; Chuikov et al., 2010; Corrigan et al., 2018).
Many ROS-related genes were altered in the Prdm16 mutant cortex (Chui et al., 2020). Chui et al. reported that changes on these genes influenced RG property in the Prdm16 mutant cortex, in which proliferation of mutant RG accelerated at E13.5 but decreased at E15.5. However, in our study we did not detect similar proliferation defects in Prdm16 mutant RG. The discrepancy could be due to the different flox lines used in these two studies: Chui et al. used a flox line that removes exon 9 of the Prdm16 gene (Cohen et al., 2014), whereas our Prdm16cGTinv (cGTinv) strain stops transcription after exon 2 (Strassman et al., 2017).
Interestingly, Hamlet, the ortholog of PRDM16, controls the temporal identity of intermediate progenitors in Drosophila neuroblast lineage (Eroglu et al., 2014), suggesting that the role of the PRDM16 proteins is evolutionarily conserved.
Transcriptional activities of PRDM16
PRDM16 can act as a repressor or an activator depending on its associated partners (Chi and Cohen, 2016). We showed that PRDM16-bound genes have a trend of derepression in the Prdm16 mutant, indicating its repressive role in the neocortex. The fact that many of the PRDM16 targets are expressed in RG suggest that repression by PRDM16 is not to fully silence genes but to maintain gene expression at the right level. In support of this, PRDM16 binding associates with open chromatin. Moreover, we do not rule out the possibility of PRDM16 being an activator, as our ChIP- and RNA-seq data identified a small subset of genes that show PRDM16 binding and downregulation in the mutant.
The PR domain of PRDM16 is essential in repressing Fezf2. Baizabal et al. showed that PRDM16 without the PR domain failed to rescue target gene derepression (Baizabal et al., 2018). The PR domain of PRDM16 was shown to be essential in suppressing MLLr1 leukemia via intrinsic H3K4 methylation activity (Zhou et al., 2016). We did not observe global changes of H3K4me3 or H3K4me2 levels in the mutant cortex by immunostaining (Fig. S10). In agreement with this result, Baizabal et al. did not detect a significant change of H3K4 methylation levels using ChIP-seq (Baizabal et al., 2018). Future studies are needed to address the mechanistic nature of how the PR domain or any other domain contributes to the function of PRDM16 in the neocortex.
Prdm16 is among the many genes deleted in human 1p36 microdeletion syndrome, a disorder that displays a wide variety of disease conditions. According to the function of PRDM16 in normal development, loss of Prdm16 might contribute to several problems including the craniofacial and cardiac defects and hydrocephalus of the syndrome (Bjork et al., 2010b; Arndt et al., 2013; Shimada et al., 2017). Our findings, along with the study by Baizabel et al., defined a mechanism by which Prdm16 loss-of-function contributes to the formation of heterotopia, a neurodevelopmental disorder that leads to severe mental retardation and seizures that were also seen in the 1p36 syndrome. More mechanistic insights of PRDM16 function in development will increase our understanding of its pathological function in diseases.
MATERIALS AND METHODS
Animals and processing
All animal procedures were approved by the Swedish agriculture board (Jordbruks Verket) with document number Dnr 11553-2017 and the MWU Institutional Animal Care and Use Committee. The Prdm16cGT and Prdm16cGTreinv mice (Strassman et al., 2017) were maintained by outcrossing with the FVB/NJ line. B6.129S2-Emx1tm1(cre)Krj/J (Emx1IREScre) (Gorski et al., 2002) were used to generate conditional gene trap knockout animals as previously described (Strassman et al., 2017).
Plasmids used for in utero electroporation include shRNA sequences targeting Flrt3 or Cdkn1c in the pll3.7-IresdsRed vector (Rubinson et al., 2003) and shPrdm16-GFP (#1) (Bjork et al., 2010a). In utero electroporation of DNA into C57J/B6L mice was performed as previously described (Saito and Nakatsuji, 2001; Saito, 2006). Briefly, timed pregnant mice with E14.5 embryos were deeply anesthetized using isoflurane and the uterine horns were exposed. We microinjected ∼1ul of plasmid DNA (3 μg/μl) spiked with Fast green (Sigma-Aldrich) into the lateral ventricles of each embryo using a glass micropipette (Drummond Scientific). Five 50 ms pulses of 40-45 mV with 950 ms intervals were delivered across the uterus with two 9 mm electrode paddles positioned on either side of the head (BTX, ECM830). As a standard procedure, the uterus was placed back in the abdominal cavity, the wound was surgically sutured and wound clips were placed. The whole procedure was performed on a warm heatpad. The health of the animals was monitored 24- and 48-h post-surgery.
The pCAGIG plasmid (Addgene plasmid #11159) was inserted with a fragment encoding a nuclear localization signal (NLS) and 3×Flag in the EcoRI site. To make pCAGIG-Prdm16-FL or pCAGIG-Prdm16-PRdeletion, the full-length open reading frame (ORF) or the truncation that lacks the coding sequence for amino acids 2-180 of Prdm16 was PCR amplified from MSCV-Prdm16 (Addgene plasmid #15504) and inserted between the EcoRI and XhoI sites in pCAGIG-NLS-Flag. The VP64 fragment was then inserted to the XhoI site of pCAGIG-Prdm16 to make pCAGIG-FL-VP64.
The plasmids used for making stable cell lines, pCDH-Prdm16 and pCDH-Prdm16-PRdeletion, were generated as follow: The Prdm16 FL ORF, the PR-deletion coding sequences or the NLS-3×Flag was digested from their pCAGIG plasmids and inserted sequentially to the pCDH-CMV-MCS-EF1-Puro plasmid (System Biosciences) between the EcoRI and NotI sites (for Prdm16-FL and Prdm16-PRdeletion) and the XbaI and EcoRI sites (for NLS-3×FLAG).
To clone shFlrt3, annealed oligos were cloned into the HpaI and XhoI sites of the lentiviral vector pLL3.7, which contains a separate CMV promoter that drives the expression of dsRed. Flrt3 shRNAs effectively suppressed co-expressed Flag-FLRT3 protein in HEK293T cells.
Immunochemistry, BrdU and EdU labelling and confocal imaging
At designed stages, embryos or pups were perfused with PBS followed by 4% paraformaldehyde. The perfused brains were dissected, fixed overnight and sectioned coronally at 60 μm using a vibratome (Leica Microsystems, VT1200S). Immunostaining was carried out according to standard protocols as previously used (Dai et al., 2013a). Details of primary and secondary antibodies used are provided in Table S4.
BrdU and EdU (5-20 μg/g of body weight) were injected into the peritoneal cavity of pregnant mice. BrdU incorporation was measured by immunostaining using an antibody against rat-BrdU (Abcam) and mouse-BrdU (Developmental Studies Hybridoma Bank; DSHB). EdU incorporation was detected with the Click-iT assay (Invitrogen) according to the manufacturer's instructions. Imaging was carried out on a Zeiss confocal microscope. ZEN (Zeiss LSM800), ImageJ (National Institutes of Health) and Photoshop (Adobe) were used for analysis and quantification.
In situ hybridization
The mouse brains at defined ages were dissected and fixed for 12 h in 4% PFA, dehydrated in 25% sucrose overnight, cryoprotected and embedded in O.C.T, and sectioned at 18 μm thickness on a Leica cryostat (CM3050s). RNA in situ hybridization was performed using digoxigenin-labeled riboprobes as previously described (He et al., 2018). Images were taken using a Leica DMLB microscope.
Quantification of in situ signal intensity was carried out using ImageJ/Fiji software. Areas at the labeled region (VZ/SVZ) and background region were selected and quantified for signal. The value was then normalized to area size.
Quantification and statistical analysis
Cell numbers were manually counted in ImageJ/Fiji cell counter. The number of marker-positive cells in the control and KO mutant at P0 was determined by counting the average number of positive cells in three 80 µm-wide columns. The number of marker-positive cells in the control and KO or cKO mutants at E15.5 was determined by counting the number of positive cells in a 100 µm-wide column from layer IV to VI. The number of marker-positive cells in the control and cKO cortex at P15 was determined by counting the number of positive cells in one 250 µm-wide column within the whole cortex in two different areas (medial and dorsal lateral). For proliferation analysis at E15.5, the number of PAX6+, TBR2+, EdU+ and KI67+ cells was determined by counting the total number of positive cells in two 100 µm-wide columns in both VZ and SVZ. For proliferation analysis at E13.5, the number of PAX6+, TBR2+, EdU+ and KI67+ cells was determined by counting the total number in the 300 µm-wide cortex. The production of daughter neurons was reflected by the cell cycle exit through measuring the ratio of EdU+KI67-cells in total EdU+ cells. The number of BrdU+CTIP2+ or BrdU−CTIP2+ cells at P5 was determined by counting the number of positive cells in a 300 µm-wide column in layer II to layer V. For in utero electroporation followed by indirect neurogenesis analysis, the number of GFP+ and TBR2+ cells was determined by counting the total number of positive cells in one 250 µm-wide column in both the VZ and SVZ. For in utero electroporation followed by migration defect analysis, the number of GFP+ cells was determined by counting the total number of positive cells in one 300 µm-wide column in both white matter and cortex. For in utero electroporation followed by progenitor proliferation analysis, the number of GFP+ and EdU+ cells was determined by counting the total number of positive cells in one 250 µm-wide column in the SVZ. All data are presented as mean±s.d., and statistical significance was determined using a two-tailed unpaired Student's t-test.
Neural stem cell culture and RT-qPCR
Control and mutant embryonic cortices were dissected and dissociated into single cell suspension and digested with Accutase (Sigma-Aldrich). Cells were maintained in proliferation media (Stemcell Technologies). Three control or three Prdm16 mutant neural stem cell cultures were grown for 2 days before RNA extraction by use of TRIzol reagent (Invitrogen). Then, 4 μg of total RNA was further cleaned with Turbo DNase (Ambion) and used in reverse-transcription with RT master mix (Thermo Fisher Scientific). To ensure the absence of genomic DNA, control qPCR was performed on a mock-reverse-transcribed RNA sample. Primer sequences are listed in Table S4.
Cell culture and luciferase assays
The neuroblastoma cell line Neuro-2A (N2A, ATCC® CCL-131™) cells were cultured in 50% of DMEM (Gibco) containing 10% fetal calf serum and 50% of Opti-MEM serum reduced medium. For luciferase assays, transfections were performed in 96-well plate using FugeneHD transfection reagent (Promega). The following DNA combinations were used: 20 ng of Fezf2, Cdkn1c, Flrt3 luciferase reporters or the pGL3 promoter vector, 100 ng of pCAGIG-Prdm16, pCAGIG-Prdm16PRdeletion, pCAGIG-Prdm16-VP64, pCAGIG-VP64 or pCAGIG. We used 2 ng of Renilla luciferase construct as an internal control. After 24 h incubation, transfected cells were lysed and luciferase activity was measured using Dual Luciferase Assays (Promega), and promoter activity was defined as the ratio between the firefly and Renilla luciferase activities.
For generating the cell lines that stably express control, PRDM16-FL or PRDM16-PRdeletion, lentiviral particles were first produced in 293T cells and then added to N2A cells for infection. The cells that stably expressed the corresponding constructs were selected and maintained in medium that contains puromycin. Two individual stable lines were generated for each of the constructs used in RT-qPCR analysis of the Fezf2 gene.
Western blots were used to examine shRNA knockdown efficiency. N2A cells were co-transfected with the constructs that express Flag-tagged FlRT3 or Flag-tagged CDKN1C and the corresponding shRNAs. After 24 h incubation, untransfected and transfected cells were harvested and lysed for protein extraction. Then the protein extract was run on 10% SDS PAGE gels. Proteins were then transferred to a PVDF membrane. The membranes were treated with blocking reagent (3% bovine serum albumin) for 1 h before incubated with the rabbit anti-RFP antibody (Abcam) overnight in a cold room. After several washes, the membranes were incubated with the secondary antibody (goat anti-rabbit immunoglobulins/HRP), washed again, and developed with the SuperSignal West Femto Maximum Sensitivity Substrate (Thermo Fisher Scientific, 34096). Images were collected using a Gel Doc XR+ Gel Documentation System (Bio-Rad) machine. The same membranes were then washed and incubated with the mouse Flag-m2 antibody (Sigma-Aldrich) overnight. On the second day, the membranes were washed, incubated with the secondary antibody (rabbit anti-mouse immunoglobulins/HRP), washed again, developed and exposed in the Gel Doc XR+ Gel Documentation System (Bio-Rad) machine. Details of the antibodies used are provided in Table S4.
In each replicate, three E13.5 control or Prdm16 KO mutant heads were pooled, fixed and lysed. ChIP was performed as previously described (Dai et al., 2013b). DNA libraries were made using the NEBNext Ultra™ II DNA Library Prep Kit and sequenced on the Illumina Hiseq2500 platform.
The replicated Prdm16 KO (×3) and control (×3) ChIP-seq samples, after the adaptor trimming by Trimmomatic, were mapped to the University of California, Santa Cruz (UCSC) Mus musculus (mm10) genome assembly using Bowtie2 with the default parameters. The uniquely mapped reads (with mapping quality ≥20) were used for further analysis. The PRDM16 peaks were called by HOMER (v4.10) (Heinz et al., 2010). The peak replicate reproducibility was estimated by IDR, using the HOMER IDR pipeline (https://github.com/karmel/homer-idr). As suggested by the Encode IDR guideline that IDR requires to initially call peaks permissively for the replicates, we used a relatively relaxed parameter ‘-F 2 -fdr 0.3 -P 0.1 -L 3 -LP 0.1’ for the true/pseudo/pooled replicates by the HOMER peak calling. The final confidence peaks were determined by an IDR<5%. The peaks that were overlapped with mm10 blacklist were also removed. For comparisons, we re-analyzed the Prdm16 control and cKO ChIP-seq public data (Baizabal et al., 2018; GSE111657) using the same HOMER IDR pipeline. All called ChIP-seq peaks are included in Table S2.
RNA-seq differential expression analysis
Cortices of control and Prdm16 KO mutant E13.5 embryos were dissected for RNA extraction using TRIzol reagent (Invitrogen). RNA quality of three biological replicates was tested by Agilent Bioanalyzer. RNA-seq libraries were made using the Illumina Truseq Total RNA library Prep Kit LT. Sequencing was performed on the Illumina Hiseq2500 platform.
After trimming the adaptor sequences using Trimmomatic, we mapped RNA-seq reads from the replicated Prdm16 wild type (×3) and mutant samples (×3) to the UCSC Mus musculus (mm10) genome assembly using HISAT2. We normalized RNA-seq by the ‘Relative Log Expression’ method implemented in the DESeq2 Bioconductor library (Love et al., 2014). Gene annotation was obtained from the iGenomes UCSC Mus musculus gene annotation. Differentially expressed mRNAs between Prdm16 mutants versus wild type were identified, and FDR (Benjamini-Hochberg) was estimated using DESeq2. For comparisons, we re-analyzed the differential expression of Prdm16 WT and cKO RNA-seq public data (Baizabal et al., 2018; GSE111660) using the same method as above. The genes with P-value ≤0.05 were considered to be differentially expressed. Differentially expressed genes are listed in Table S1.
The ATAC-seq libraries were made according to the published method (Buenrostro et al., 2013) and using the Illumina Nextera DNA library kit. In brief, cortices were dissected from three control and three Prdm16 KO E13.5 brains. Tn5 enzyme reaction was performed at 37 C for 30 min, followed by DNA purification. Eleven cycles of PCR amplification were performed using barcoded adaptors and primers on a purified DNA template. Libraries were purified and pooled before sequencing with the Illumina Next-seq platform. The replicated Prdm16 KO (×3) and control (×3) ATAC-seq samples, after the adaptor trimming by Trimmomatic, were mapped to the UCSC Mus musculus (mm10) genome assembly using Bowtie2 with the default parameters. The high quality and uniquely mapped reads (with mapping quality ≥20) were used for further analysis. ATAC-seq differential expression analysis between Prdm16 mutants and wild types on the Prdm16-bound ChIP-seq peaks was performed using the Limma R package. The ATAC-seq peak calling was performed by HOMER using the ‘broad peak’ option with parameters ‘-region -size 1000 -minDist 2500’, separately for the mutant and wild type. To compare active enhancers between E13.5 and E15.5, we further re-analyzed the publicly available histone mark H3K27ac and PRDM16 ChIP-seq data at E15.5. We called the peaks against input using ‘narrow peak’ option in HOMER with the default parameters.
Gene set enrichment testing
To test whether a set of genes is significantly changed amongst the differentially expressed genes from Prdm16 wild type and mutant RNA-seq data, we used gene set testing function ‘camera’ and ‘mroast’ in the Limma R package (Ritchie et al., 2015). We used ‘camera’, a ranking-based gene set test accounting for inter-gene correlation, to test whether the layer markers were significantly changed as a set. We used ‘mroast’ (number of rotations=1000), a self-contained gene set test, to test whether the majority of the genes amongst PRDM16 targets were significantly up- or downregulated. We also used ‘mroast’ to test which GO terms and Reactome pathways were significantly up- or downregulated in Prdm16 mutant versus wild type.
To gain insights into cell types of the Prdm16 targets, we reanalyzed the murine cortical time-series scRNA-seq data (Yuzwa et al., 2017). We employed the Bioconductor scRNA-seq analysis workflow for droplet-based protocols (Lun et al., 2016). (1) The cortical cells (the cells expressing Emx) were selected for the analysis. The low quality cells were first removed if they were 3 median absolute deviation (MAD) lower than the median library size or if they were 3 MAD lower the median gene expression or if they were 4 MAD higher than the median mitochondrial reads. (2) We used the deconvolution approach, a method to handle high zeros in scRNA-seq, to compute size factors for cells for normalization. (3) The cells were constructed into graphs by creating a shared nearest neighbor graph and clustered by Walktrap algorithm. (4) We manually assigned cell types to the identified clusters in each stage using the known neuron and layer markers. For each cell, the sum of (log2) reads for all Prdm16 targets was calculated. The number of cells with >200 sum of reads divided by the total number of cells in each cluster were calculated to give the percentage of cells that Prdm16 targets expressed in Fig. 6.
To identify differentially expressed genes between E13.5 and E15.5, RG cells were extracted from E13.5 and E15.5 RG clusters and differential expression analysis was performed using the edgeRQLF R package. The genes with FDR≤0.2 and FC>1.4-fold were considered to be significantly differentially expressed between E13.5 and E15.5, and are listed in Table S3.
We thank the animal experimental core facility and the imaging facility of Stockholm University and the National Genomic Institute of Scilife Laboratories, Sweden, for their service and support.
Conceptualization: Q.D.; Methodology: L.H., J.J., W.H., B.C.B., J.W., Q.D.; Software: J.W.; Validation: L.H., J.J., W.H.; Formal analysis: L.H., J.W., Q.D.; Investigation: L.H., J.W.; Resources: L.H., J.J., B.C.B.; Data curation: L.H., J.W.; Writing - original draft: Q.D.; Writing - review & editing: L.H., B.C.B., J.W., Q.D.; Visualization: L.H.; Supervision: B.C.B., Q.D.; Project administration: Q.D.; Funding acquisition: Q.D.
W.H. was supported by the visiting scientist fellowship from the China Scholarship Council. J.W. was supported by the Australian Research Council (ARC) Future Fellowship (FT60100143). Work from the B.C.B. lab was supported by the Midwestern University Office of Research and Sponsored Programs for intramural funding. The project was supported by the Young Investigator grant from Swedish Research Council (Vetenskapsrådet, 2014-5584) and the research grant from Swedish Cancerfonden (CAN 2017/529) to Q.D. Deposited in PMC for immediate release.
All deep sequencing data produced in the study have been deposited in GEO under accession number GSE162854.
Peer review history
The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.194670.reviewer-comments.pdf
The authors declare no competing or financial interests.