Single-Cell Analysis of Foxp1-Driven Mechanisms Essential for Striatal Development.

The striatum is a critical forebrain structure integrating cognitive, sensory, and motor information from diverse brain regions into meaningful behavioral output. However, the transcriptional mechanisms underlying striatal development at single-cell resolution remain unknown. Using single-cell RNA sequencing (RNA-seq), we examine the cellular diversity of the early postnatal striatum and show that Foxp1, a transcription factor strongly linked to autism and intellectual disability, regulates the cellular composition, neurochemical architecture, and connectivity of the striatum in a cell-type-dependent fashion. We also identify Foxp1-regulated target genes within distinct cell types and connect these molecular changes to functional and behavioral deficits relevant to phenotypes described in patients with FOXP1 loss-of-function mutations. Using this approach, we could also examine the non-cell-autonomous effects produced by disrupting one cell type and the molecular compensation that occurs in other populations. These data reveal the cell-type-specific transcriptional mechanisms regulated by Foxp1 that underlie distinct features of striatal circuitry.

INTRODUCTION

The striatum is the major input nucleus of the basal ganglia and receives dense glutamatergic inputs from the cortex and thalamus, as well as dopaminergic innervations from the substantia nigra (SN) and other neuromodulatory circuits. The principal neurons that receive and integrate this information within the striatum are GABAergic spiny projection neurons (SPNs) (Gerfen and Surmeier, 2011). Proper function of striatal circuitry is essential for coordinated motor control, action selection, and rewardbased behaviors (Cui et al., 2013; Tecuapetla et al., 2016). Dysfunction of this system is implicated across many neurological disorders, including Huntington’s disease, Parkinson’s disease, autism spectrum disorder (ASD), and obsessive-compulsive disorder (Crittenden and Graybiel, 2011; Fuccillo, 2016).

Striatal organization has two prominent features: the division of the striatum into distinct neurochemical zones, the striosome and matrix compartments, and the division of SPNs into the direct or indirect projection pathways. Striosome and matrix compartments are enriched for distinct neuropeptides and contribute differentially to striatal connectivity and behavior (Crittenden and Graybiel, 2011; Crittenden et al., 2016; Friedman et al., 2015; Smith et al., 2016). Recent evidence suggests that striosome-matrix compartmentalization is the initial organizational plan during striatal development, with distinct intermediate progenitor pools in the lateral ganglionic eminence (LGE) giving rise first to striosome SPNs and then to matrix SPNs (Kelly et al., 2018). These progenitor pools then generate either direct or indirect pathway SPNs, which populate both compartments (Kelly et al., 2018). Direct pathway SPNs (dSPNs) express dopamine receptor 1 (D1) and project to the globus pallidus internal (GPi) and SN. Indirect pathway SPNs (iSPNs) express dopamine receptor 2 (D2) and project to the globus pallidus external (GPe) (Gerfen and Surmeier, 2011). Ultimately, these pathways work to bidirectionally modulate excitatory inputs back onto the cortex (Gerfen and Surmeier, 2011).

Mature dSPNs and iSPNs have distinct molecular profiles based on expression profiling studies (Gokce et al., 2016; Heiman et al., 2008; Lobo et al., 2006; Saunders et al., 2018), and several transcription factors and chromatin regulators have been identified for both pan-SPN and dSPN/iSPN sub-specification (Anderson et al., 1997; Arlotta et al., 2008; Corbin et al., 2000; Ehrman et al., 2013; Garel et al., 1999; Kim et al., 2008; Lobo et al., 2008; Long et al., 2009; Lu et al., 2014; Martín-Ibáñez et al., 2017; Maze et al., 2014; Waclaw et al., 2017; Xu et al., 2018). However, the molecular mechanisms governing both SPN specification and striosome-matrix organization remain incomplete.

Forkhead-box protein 1 (Foxp1) is a transcription factor with enriched expression in the striatum compared to the rest of the brain (Heiman et al., 2008). Expression of Foxp1 begins in the marginal zone of LGE between embryonic day 12 (E12) and E13 and is maintained throughout striatal development in both SPN subtypes (Ferland et al., 2003; Precious et al., 2016). Loss-of-function FOXP1 variants are strongly linked to ASD and intellectual disability in humans (Fong et al., 2018; Meerschaut et al., 2017; Siper et al., 2017). We previously found that Foxp1 was critical for maintaining the intrinsic excitability of iSPNs within the striatum of Foxp1 heterozygous mice (Araujo et al., 2015), and two studies have suggested a general role for Foxp1 in striatal development (Bacon et al., 2015; Precious et al., 2016). However, no study has examined the contribution of Foxp1 to striatal development with cell-type specificity or the downstream targets of Foxp1 at single-cell resolution. Given our previous results, we hypothesized that Foxp1 differentially regulates the development, molecular composition, and function of distinct SPN subpopulations.

To ascertain the cell-type-specific role of Foxp1, we generated mice with deletion of Foxp1 from dSPNs, iSPNs, or both populations and used a combination of single-cell RNA sequencing (scRNA-seq), serial two-photon tomography, and behavioral assays to delineate the contribution of Foxp1 to striatal development and function. We show that Foxp1 is crucial for maintaining the cellular composition of the striatum, especially iSPN specification, and proper formation of the striosome-matrix compartments. We uncover downstream targets regulated by Foxp1 within iSPNs and dSPNs and connect these molecular findings to cell-type-specific deficits in motor and limbic system-associated behaviors, including motor learning, ultrasonic vocalizations (USVs), and fear conditioning. We overcome a limitation of previous studies by identifying the non-cell-autonomous gene expression changes and functional consequences when distinct SPN subtypes are manipulated. These findings provide an important molecular window into postnatal striatal development, which has not been characterized at single-cell resolution, and further our understanding of striatal circuits mediating ASD-relevant behavioral phenotypes.

RESULTS

Early Postnatal scRNA-Seq of Striatal Cells across Foxp1 cKO Mice

To examine the contribution of Foxp1 to striatal development in a cell-type-specific manner, we generated Foxp1 conditional knockout (cKO) mice using bacterial artificial chromosome/ clone (BAC)-transgenic mice driving Cre expression under the D1- or D2-receptor promoters (Gong et al., 2007) crossed to Foxp1 mice (Araujo et al., 2017; Feng et al., 2010; Usui etal.,2017)(Figure1A).Fourgenotypeswereusedfordownstream analyses: Drd1-Cre; Foxp1 (Foxp1, deletion of Foxp1 in dSPNs), Drd2-Cre ; Foxp1 (Foxp1, deletion of Foxp1 in iSPNs), Drd1-Cre ; Drd2-Cre ; Foxp1 (Foxp1, deletion of Foxp1 in both d/iSPNs), and Foxp1 (Foxp1). We confirmed that Foxp1 was reduced at both the transcript and protein levels within the striatum at postnatal day 7 (P7) and P56 (Figures 1B–1D). By crossing the Drd1- and Drd2-Cre lines to R26-stop-eYFP mice (Figure S1A), we found that both Cre-lines turn on and co-localize with Foxp1 within a 24-h window of each other between E14 and E15 (Figures S1B–S1G). While there is expression of D1 and D2 receptors outside of the striatum, there are few regions of notable overlap with Foxp1+ cells outside of lower-layer cortical neurons that highly express Drd1 (Figures S1H and S1I).

Figure 1.

Early Postnatal scRNA-Seq of Striatal Cells across Foxp1 cKOs

(A) Schematic of the scRNA-seq experiment using striatal tissue from P9 mice (n = 4/genotype) with cell-type-specific conditional deletion of Foxp1 within the dopamine receptor 1 (Foxp1), dopamine receptor 2 (Foxp1), or both (Foxp1) cell types.

(B–D) Foxp1 is reduced in the striatum via immunohistochemistry (P7 and P56) (B) and quantitative RT-PCR (P7) (C) within each cKO line, with near-complete reduction in Foxp1 striatal tissue via immunoblot (P7) (D) (scale bar, 100 μm).

(E) Non-linear dimensionality reduction with UMAP of all 62,778 post-filtered cells combined across genotype and used for downstream analyses. Cell-type annotation is overlaid to identify the major cell type represented by each cluster (43 total clusters).

(F) Stacked bar plots of the contribution of cells from each genotype to each cluster, with arrows indicating genotype-driven changes within SPN (blue arrow) orneurogenic progenitor (purple arrow) clusters.

(G) Pie charts using colors from (E) show the striatal cell-type composition as a percentage of total cells within each genotype.

See also Figures S1 and S2 and Tables S1 and S2.

Using 10X Genomics Chromium technology (Zheng et al., 2017), we profiled the transcriptome of 62,778 striatal cells across control and the three Foxp1 cKO mouse lines at P9 (n = 4/genotype; 16 samples total) (Figure 1E). This time point is an important and understudied period of striatal development before excitatory synaptic density onto SPNs markedly increases and where perturbations of cortical-striatal activity can have long-lasting effects on SPN spine density and circuit activity (Kozorovitskiy et al., 2012; Peixoto et al., 2016). We detected 5,587 unique molecular identifiers (UMIs) (median = 3,837) and 1,794 genes (median = 1,532) per cell across all genotypes (Figure S2A). All cells were combined across genotype and filtered for downstream clustering, resulting in 43 clusters driven primarily by cell type (Figures 1E and 1F; Table S1). For unbiased characterization of striatal cell types, we used a previously annotated adult striatal single-cell dataset (Saunders et al., 2018) to assign cell types to each cluster using two separate methods, a previously published expression weighted cell-type enrichment (EWCE) analysis (Skene and Grant, 2016) and an in-house correlation analysis (see STAR Methods). We confirmed the cell-type annotation of our dataset by examining the expression of known marker genes for each major cell type and found similar cell-type composition variability between replicates within each genotype (Figures S2B–S2D; Table S1). The principal cell types found within the early postnatal striatum were SPNs, neurogenic progenitor cells, astrocytes, and oligodendrocyte precursor cells (OPCs) (Figure 1G; Table S2). Endothelial, microglia, ependymal, interneurons, and mural cells made up a smaller percentage of total cells within the postnatal striatum (Figure 1G). We found clusters expressing proliferation markers (Mki67; Figure S2E), progenitor markers (Ascl1 and Dlx2; Figures S2E and S2F), and SPN-specification markers (Sp9, Ppp1r1b, Drd1, and Drd2; Figure S2F).

The cell type with the largest number of subclusters were SPNs, with 13 distinct clusters (Figure 1E). SPNs and neurogenic progenitors made up 52% of the total cell population (Figure 1G), and genotype-specific variations were observed primarily within these clusters, where Foxp1 is selectively deleted (arrows in Figure 1F; Table S2). To more directly compare the composition of striatal cell types across genotypes and better control for variations in total cells sequenced between genotypes, we down-sampled the dataset to yield equal cell numbers across genotypes and reclustered the resultant cells separately. We found analogous resultsin the percentage of cell types from down-sampling experiments compared to the full dataset (Figure S1G).

Diversity of Early Postnatal Striatal Projection Neurons

To further characterize early postnatal SPN subtypes and the effects of Foxp1 deletion, we next isolated all clusters identified as neuronal from the annotation analyses (see Figure 1; STAR Methods) and reclustered them separately (18,073 cells total and 24 clusters) (Figure 2A). We also labeled each cell by genotype (Figures 2B–2E) and examined SPN-subtype identity. Three interneuron clusters (clusters 14, 15, and 20) were identified by the interneuron marker Nkx2–1 (Figures 2A and S3A; Table S1). We could clearly distinguish dSPN clusters (clusters 0, 1, 3, 4, 5, and 9) and iSPN clusters (clusters 2, 8, 10, and 16) using canonical markers (Drd1 and Tac1 for dSPNs and Drd2 and Penk for iSPNs) (Figures 2F, S3B, and S3C; Table S1). Pairwise comparisons between the major dSPN and iSPN clusters confirmed enrichment of known genes within each population (Figure S3D). One small cluster (cluster 19) co-expressed both Drd1 and Drd2 receptors, termed “ddSPNs” (Figure 2F). SPNs expressing both Drd1 and Drd2 receptors were also scattered throughout other clusters and comprised ~1% of the total SPN population (Figure S3E; Table S2). We identified a recently described SPN subpopulation termed “eccentric” SPNs (eSPNs) (Saunders et al., 2018) within cluster 7 that distinctly expressed markers, such as Casz1 and Otof (Figures 2F and S3F; Table S1). We found two clusters (clusters 6 and 23) that were enriched for the neurogenic transcription factors Sox4 (Figures 2F and S3G) and Sox11 (Table S1). Sox4 and Sox11 function during the terminal steps of neurogenesis to promote neuronal maturation (Bergsland et al., 2006; Chen et al., 2015); therefore, we termed these clusters “immature” SPNs (imSPNs). We confirmed the presence of Sox4+ cells within and near the subventricular zone of the lateral ventricle and populating zones in P7 ventral striatum (Figures S3G–S3H). Additionally, several clusters enriched for d/iSPN markers also have high expression of Sox4 (dSPN clusters 9, 11, 13, and 17 and iSPN cluster 16), indicating these may be less mature SPNs (Figures 2F and S3G). Two clusters (clusters 12 and 18) were composed primarily of cells from Foxp1 and Foxp1 and could not be classified into distinct SPN subclusters, termed “unknown” (Figures 2A–2F). Foxp2, another Foxp transcription factor with high sequence similarity to Foxp1 (Shu et al., 2001), is an SPN marker with enriched expression in dSPNs (Figures S3D and S3I) (Fong et al., 2018; Vernes et al., 2011). Within our dataset, Foxp2 is highly expressed within all dSPN clusters and one iSPN cluster (cluster 8). Surprisingly, the highest expression of Foxp2 is found within eSPN cluster 7 and imSPN cluster 6, where notably Foxp1 is not highly expressed (Figures 2F and S3J). Foxp2 expression is also maintained within adult eSPNs (Saunders et al., 2018). We confirmed that Foxp2 is expressed in cells other than D1- and D2-expressing cells at P9 using D1-tdTomato and D2-eGFP reporter mice (Figure S2K).

Figure 2.

Foxp1 Specifies Distinct SPN Subpopulations

(A) UMAP plot showing each neuronal subcluster by color with overlay showing neuronal subpopulation identity.

(B–E) UMAP plots of cells from (A) color-coded to identify each cell by genotype.

(F) Violin plots of the normalized UMI expression of markers of SPN subpopulations: dSPNs (Drd1, Tac1, and Foxp2), iSPNs (Drd2 and Penk), ddSPNs (Drd2, Drd1, and Tac1), eSPNs (Casz1), and imSPNs (Sox4).

(G) Stacked bar plots of the contribution of cells from each genotype to each cluster. Cluster numbers in red indicate dSPNs, blue indicate iSPNs, and italicizednumbers indicate Foxp1 cKO driven SPN clusters.

(H) Pie charts showing altered composition of SPN subtypes within Foxp1 cKO mice (using colors from A).

(I and J) Foxp1 cKO mice were crossed to D2 reporter lines to label dopamine receptor 2 (D2) iSPNs in green (coronal section, scale bar, 500 μm) (I). Foxp1 and Foxp1 mice had significantly fewer iSPNs compared to Foxp1 mice at P7, while Foxp1 mice had significantly more iSPNs compared to Foxp1 animals. Data are represented as a boxplot (J); n = 3–6 mice/genotype; ****p <0.0001 and *** p <0.005, one-way ANOVA with Tukey’s multiple comparisons test.

See also Figure S3 and Tables S1 and S2.

Foxp1 Regulates SPN Subtype Composition

We next asked whether Foxp1 regulates the development of specific SPN subtypes. We observed that dSPNs and/or iSPNs with deletion of Foxp1 largely clustered separately from control SPNs (Figures 2B–2E and 2G). By examining the percentage of SPN subtypes across genotypes, we found that control samples have nearly double the number of dSPNs relative to iSPNs (61% dSPNs and 31% iSPNs), with imSPNs contributing ~4% of the total SPN population and both eSPNs and ddSPNs contributing ~2% (Figure 2H). This percentage of dSPNs to iSPNs at P9 is similar to those seen at P14 using reporter mice (Thibault et al., 2013). The percentage of SPN subtypes varied across Foxp1 cKO samples (Figure 2H). Notably, the number of eSPNs increased 2- to 4-fold across Foxp1 cKO samples and, within Foxp1 and Foxp1 samples, the percentage of iSPNs was reduced by two-thirds compared to controls (Figure 2H).

To independently confirm the reduction of iSPNs in Foxp1 and Foxp1 samples, we crossed all Foxp1 cKO mice to D2-eGFP reporter animals (D2-eGFP; Foxp1) to label iSPNs (Figure 2I). Within Foxp1 mice, we again found a significant two-thirds reduction of iSPNs (D2-eGFP+ cells) as seen in the scRNA-seq data (Figures 2I and 2J). Compared to Foxp1, Foxp1DD mice also showed significantly reduced iSPNs, but they also showed increased iSPNs compared to Foxp1 (Figures 2I and 2J). The remaining iSPNs in the Foxp1 mice were not the product of D2-Cre inefficiency, as these cells did not express Foxp1 (Figure S3L). Only seven iSPNs within the single-cell Foxp1 data did not express Foxp1 (0.2% of total iSPNs) (Figure S3E); therefore, we would not expect the remaining iSPNs in Foxp1D2 mice to be a naturally occurring Foxp1-negative population. Taken together, these results indicate that Foxp1 is required for the development of a distinct iSPN subpopulation.

Deletion of Foxp1 Disrupts iSPN Localization to Striosomes and Striosomal Area

To determine whether the remaining subpopulation of iSPNs within Foxp1 or Foxp1 mice localized to the striosome or matrix compartment, we stained for the canonical striosome marker MOR (Oprm1) in Foxp1 cKO mice crossed to D2-eGFP reporter mice (Figures 3A–3H). We found that few remaining iSPNs within Foxp1 and Foxp1 mice localized within the striosome compartment compared to control and Foxp1 mice (Figure 3I). Instead, the remaining iSPNs clustered around the border of striosomes, with some iSPNs scattered throughout the matrix (Figures 3A–3H). We also found reduced striosomal area across all Foxp1 cKO animals at P7 (Figure 3J), and fewer striosome “patches” were observed within Foxp1DD mice (Figure 3K). This reduction in striosomal area was also found using a different striosomal marker in early postnatal development, PPP1R1B (also known as DARPP32) across Foxp1 cKO mice (Figures S4A and S4B).

Figure 3.

Foxp1 Regulates Striosome-Matrix Organization

(A–H) Immunohistochemistry for mu-opiod receptor (MOR) in P7 striatal sections from Foxp1 (A, E), Foxp1 (B, F), Foxp1 (C, G) and Foxp1 (D, H) mice crossed to D2-eGFP reporter mice to label D2+ (EGFP) cells (scale bars represent 500 μm in A–D and 100 μm in E–H). (E–H) Zoomed in images of white-dashed squares in (A)-(D).

(I) Quantification of the number of D2+ cells in striosomes (MOR+) relative to the total number of D2+ cells across control and Foxp1 cKO mice. Data are represented as mean ± SEM; n = 5–9 slices/2 mice/genotype; ***p < 0.0001, one-way ANOVA with Dunnett’s multiple comparisons test.

(J and K) The striosome compartment was significantly reduced across all Foxp1 cKO mice as a percent of total striatal area (measuring only dorsal striosomes) (J), and the number of striosome “patches” was significantly reduced in Foxp1 animals (K). Data are represented as mean ± SEM; n = 9 slices/3 mice/genotype; *p < 0.05, p** < 0.005, and ***p <0.0001, one-way ANOVA with Dunnett’s multiple comparisons test.

(L) Scatterplots showing the percent expression of enriched transcripts between Foxp1 iSPN cluster 2 versus cluster 8.

(M) Scatterplots showing pairwise comparison of percent expression of enriched transcripts between Foxp1 dSPN clusters (cluster 0 versus cluster 5, cluster 0 versus cluster 9, and cluster 5 versus cluster 9).

(N) Scatterplots showing pairwise comparison of percent expression of enriched transcripts between Foxp1 dSPN clusters (cluster 1 versus cluster 3, cluster 1 versus cluster 4, cluster 3 versus cluster 4). Striosome and matrix markers are indicated within each scatterplot (p.adj < 0.05, percent expression > 0.2).

See also Figure S4 and Table S3.

Using a pairwise differential gene expression analysis between clusters in our scRNA-seq data, we found subclusters within mature dSPNs and iSPNs enriched for known markers of striosome or matrix compartments (Figure S4C). We found two Foxp1 iSPN clusters enriched for known markers of the matrix (cluster 2, Penk, Chrm3, and Epha4) or striosome (cluster 8, Nnat, Lypd1, and Foxp2) compartments (Figure 3L). Remaining iSPNs within Foxp1 and Foxp1 samples clustered distinctly into one cluster (cluster 10) enriched for matrix markers Penk and Chrm3, confirming our immunohistochemistry (IHC) data showing significant loss of striosomal iSPNs with deletion of Foxp1 (Table S1). In dSPN clusters, we found three Foxp1 dSPN clusters (clusters 0, 5, and 9) (Figure 3M). Cluster 0 was enriched for matrix markers (Ebf1, Epha4, and Mef2c), and cluster 5 was enriched for striosome markers (Oprm1, Isl1, Pdyn, Lypd1, Tac1, and Nnat) (Figure 3M; Table S3). Cluster 9 had both matrix (Epha4, Brinp3) and striosome enriched markers (Sepw1 and Ppp1r1b) (Figure 3M; Table S3). Within Foxp1 dSPN clusters (clusters 1, 3, and 4), there is a loss of distinct striosome-matrix molecular enrichment (Figures 3N and S4C). These data indicate that Foxp1 plays an important role within both dSPNs and iSPNs to maintain proper striosome-matrix architecture and cellular composition.

Cell-Type-Specific Foxp1 Regulated Targets

To better understand the molecular mechanisms regulated by Foxp1, we performed a cell-type-specific “pseudobulk” differential gene expression analysis (see STAR Methods) of the scRNA-seq data across genotypes. We identified differentially expressed genes (DEGs) regulated by Foxp1 within dSPNs or iSPNs, both cell-autonomously and non-cell-autonomously (Figures 4A and 4B; Table S4). Cell-autonomous DEGs are found in Cre active cells (dSPNs in Foxp1 samples or iSPNs in Foxp1 samples) and non-cell-autonomous DEGs are found in Cre inactive cells (iSPNs in Foxp1 samples or dSPNs in Foxp1 samples). We observed more total iSPN-DEGs (647) compared to dSPNs-DEGs (285) across genotypes (Figures 4A and 4B). There were more cell-autonomous changes than non-cell-autonomous changes within both dSPNs and iSPNs of Foxp1 and Foxp1 samples, and no differences in the ratio of cell-autonomous to non-cell autonomous DEGs within dSPNs or iSPNs were observed (Figure 4C). However, significantly more iSPN-DEGs were shared between Foxp1 and Foxp1 samples (211 DEGs) compared to dSPN-DEGs shared between Foxp1 and Foxp1 samples (47 DEGs) (Figure 4D). The DEGs specific to Foxp1DD samples were termed “interaction” DEGs. We found significantly more interaction DEGs in dSPNs, suggesting that iSPN dysfunction exerts more transcriptional changes upon dSPNs with loss of Foxp1 than vice versa (Figure 4D).

Figure 4.

Foxp1 Regulates Cell-Type-Specific Molecular Pathways

(A and B) SPN cell-type-specific differential gene expression between genotypes. Upset plot showing the overlap of upregulated or downregulated DEGs across genotypes within iSPNs (A) or dSPNs (B). Genes shown within boxes are color-coded by categories indicated.

(C) No significant difference between the number of DEGs within iSPNs and dSPNs that are cell-autonomous versus non-cell-autonomous (p = 0.0975, two-sidedFisher’s exact test).

(D) There is a significant difference in the number of DEGs within Foxp1 mice that overlap with Foxp1 or Foxp1 DEGs to specific Foxp1 DEGs (interaction DEGs) (p < 0.0001, two-sided Fisher’s exact test).

(E and F) Enrichment of ASD-risk genes SFARI score 1–4 with upregulated or downregulated iSPN-DEGs (E) blue) or dSPN-DEGs (F), red) across Foxp1 cKO samples using a hypergeometric overlap test (8,000 genes used as background).

(G and H) Enrichment of upregulated or downregulated iSPN-DEGs (G) or dSPN-DEGs (H) across Foxp1 cKO samples in distinct SPN subtypes (top 50 most enriched genes/cluster) using a hypergeometric overlap test (8,000 genes used as background).

See also Figure S5 and Tables S4 and S5.

The difference in total number of DEGs between iSPNs and dSPNs could be due to transcriptional compensation by Foxp2 in dSPNs. Foxp2 is enriched in dSPNs relative to iSPNs (Figure S2C), and we previously found that Foxp1 and Foxp2 have shared striatal targets (Araujo et al., 2015). Interestingly, Foxp2 is increased in iSPNs with loss of Foxp1, suggesting that Foxp1 may function to repress Foxp2 within distinct iSPN subtypes (Figure 4A; Table S4). Six3 (Six homeobox 3), a transcription factor crucial for iSPN specification (Xu et al., 2018), is also upregulated within the remaining iSPNs of Foxp1 and Foxp1 mice (Table S4). We previously found that SIX3 was a direct target of FOXP1 in human neural progenitors (Araujo et al., 2015). Therefore, upregulation of both Foxp2 and Six3 in iSPNs may play a role in the specification of the remaining iSPNs within Foxp1 and Foxp1 mice.

Gene Ontology (GO) analysis of the shared iSPN-DEGs within Foxp1 and Foxp1 supports a role for Foxp1 in axon guidance, neurogenesis, and neuronal differentiation of iSPNs (Table S5). Shared upregulated dSPN-DEGs within Foxp1 and Foxp1 suggest altered synaptic and voltage-gated mechanisms (Table S5). We confirmed changes in cell-type-specific gene expression via IHC for a subset of top DEGs (Pde1a, Calb1, and Ppp1r1b) using dual-reporter mice labeling dSPNs with tdTomato (Drd1-tdTomato; Foxp1) and iSPNs with eGFP (Drd2-eGFP; Foxp1) crossed to Foxp1 cKO strains (Figure S5). Pde1a, a gene encoding a calmodulin/ Ca2+-activated phosphodiesterase, was upregulated in both SPN subtypes within all Foxp1 cKO samples in a cell-autonomous and non-cell-autonomous manner (Figures 4A, 4B, S5A, and S5D–S5F). Previous in vitro work found that loss of Foxp1 reduced the expression of PPP1R1B, a critical phosphatase in the dopamine signaling cascade (Precious et al., 2016). We show this decrease in PPP1R1B is specific to iSPNs in vivo (Figures S5B and S5D–S5F). We also confirmed the increase of calbindin 1 (Calb1) selectively in dSPNs with deletion of Foxp1 (Figures S5C–S5F).

Given our previous finding that striatal targets of Foxp1 overlapped significantly with ASD-associated genes (Araujo et al., 2015), we examined the cell-type specificity of this overlap (Figures 4E and 4F). Using the SFARI ASD gene list, we found a significant overlap with high-confidence ASD-risk genes (SFARI gene score 1–4) with iSPNs-DEGs with cell-autonomous deletion of Foxp1. These genes included three members of the contactin family of axon-associated cel-adhesion molecules: Cntn4, Cntn5, and Cntn6 (Figures 4A and 4E). There was no significant overlap with ASD-risk genes and cell-autonomous DEGs in dSPNs (Figure 4F). Surprisingly, we found a significant overlap with upregulated, non-cell-autonomous iSPN-DEGs (four DEGs total) within Foxp1 samples that were ASD genes (Kirrel3 and Nlgn1) (Figure 4E). Both iSPN- and dSPN-DEGs within Foxp1 samples overlapped with ASD-risk genes (Figures 4E and 4F). These data demonstrate that cell-type-specific deletion of Foxp1 specifically within iSPNs modulates ASD-associated molecular pathways both cell autonomously and non-cell autonomously.

Two ASD-risk genes that were upregulated with deletion of Foxp1 in dSPNs were Cntnap2 (contactin-associated protein like 2) and Dpp10 (dipeptidyl peptidase-like 10) (Figure 4B; Table S4). Cntnap2 is a known repressed downstream target of both Foxp1 and Foxp2 (O’Roak et al., 2011; Vernes et al., 2008), and we previously found upregulation of Dpp10 within Foxp1 striatal tissue using bulk RNA-seq (Araujo et al., 2015). Here, using scRNA-seq, we show this regulation is specific to dSPNs.

Upregulation of eSPN Molecular Markers with Deletion of Foxp1

To determine whether deletion of Foxp1 within SPNs altered cell identity, we overlapped the top 50 enriched gene markers of distinct SPN subpopulations (eSPNs, imSPNs, and matrix and striosome dSPNs and iSPNs) (Table S1) with upregulated or downregulated iSPN-DEGs (Figure 4G) or dSPN-DEGs (Figure 4H) found within each Foxp1 cKO group. The upregulated DEGs in both iSPNs and dSPNs with cell-autonomous deletion of Foxp1 were significantly enriched for molecular markers of eSPNs (Figures 4G and 4H). Upregulated iSPN-DEGs were specifically enriched for the top four enriched eSPNs markers (Adarb2, Ntng1, Asic2, and Foxp2 (Figure 4A). We also performed a pseudobulk DEG analysis removing eSPNs and still found significant enrichment of eSPN markers upregulated in both iSPNs and dSPNs (Figures S5G and S5H). iSPN and dSPN subtype enriched genes significantly overlapped with downregulated DEGs in both Foxp1 and Foxp1 samples (Figures 4A, 4G, and 4H). Taken together, these results indicate that Foxp1 is important for maintaining the molecular identity of dSPNs and iSPNs within both matrix and striosome compartments and repressing the molecular identity of a recently described potential SPN subtype, eSPNs.

Altered Direct Pathway GPi Collaterals and Indirect Pathway GPe Projections in Foxp1 Mice

Many DEGs regulated by Foxp1 within SPNs are involved in axonogenesis and neuron projection (Tables S4 and S5). We therefore examined SPN projection patterns impacted by cell-typespecific deletion of Foxp1 in adult mice using serial two-photon tomography combined with a machine-learning-based quantification algorithm (Ragan et al., 2012; Sommer et al., 2011). We crossed Foxp1, Foxp1, and Foxp1 mice to D1-tdTomato and/or D2eGFP reporter mice (described above) to visualize projection patterns of both the direct (dSPN) and indirect (iSPN) pathway, respectively (Figures 5A–5C). We found a significant reduction of iSPN terminals onto the GPe in Foxp1 mice, which was not unexpected given the reduced number of iSPNs (Figures 5B and 5D). iSPN terminals onto the GPe were unaltered in Foxp1 mice (Figures 5C and 5D). Moreover, there were no changes in dSPN projection patterns in Foxp1 mice (Figures 5C and 5E); however, Foxp1 mice had significant deficits in dSPN collaterals onto the GPi, supporting a non-cell-autonomous role for Foxp1 in iSPNs (Figures 5B and 5E). We also quantified total striatal area across genotypes and found a significant decrease in striatal area in Foxp1 mice, while no changes were found in Foxp1 animals (Figure 5F). These findings indicate that Foxp1 regulates both iSPN and dSPN projection patterning through its role in iSPNs (Figure 5G).

Figure 5.

Deletion of Foxp1 in iSPNs Alters Projection Patterns of Both dSPNs and iSPNs

(A–C) Representative Tissuecyte 1000 coronal sections (top panels) and 3D image of probability map projected onto reference brain (gray) (bottom panels) showing the projections of dSPNs and iSPNs using D1 and D2 reporter mice, respectively, crossed to Foxp1 (A), Foxp1 (B), or Foxp1(C).

(D) Quantification of the normalized probability maps of iSPN (EGFP) projections within Foxp1, Foxp1, and Foxp1 mice showing reduced GPe projections from iSPNs within Foxp1 mice. No significant changes were seen in projection patterns onto the SNc or SNr. Data are represented as mean ± SEM; n = 3–4 mice/genotype; ***p <0.0001, two-way ANOVA with Dunnett’s multiple comparisons test.

(E) Quantification of the normalized probability maps of dSPN (tdTomato) projections within Foxp1, Foxp1, and Foxp1 mice showing reduced GPi projections from dSPNs within Foxp1 mice. Data are represented as mean ± SEM; n = 2–4 mice/genotype; **p <0.01, two-way ANOVA with Dunnett’s multiple comparisons test.

(F) Striatal area quantification of four serial slices from anterior to posterior at 400-μm increments within Foxp1, Foxp1, and Foxp1 adult mice. Data are represented as mean ± SEM; n = 3–4 mice/genotype; ***p <0.001, one-way ANOVA with Dunnett’s multiple comparisons test.

(G) Schematic of cell-autonomous and non-cell-autonomous projection deficits found to the GPe and GPi in Foxp1 animals.

(H) Overlap of dSPN-DEGs within Foxp1 or Foxp1 cells. Foxp1 dSPN-DEGs that are involved in neuron projection are shown, with ASD-risk genes highlighted in purple.

GPe, globus pallidus external; GPi, globus pallidus internal; SNc, substantia nigra pars compacta; SNr, substantia nigra pars reticulata; STR, striatum.

Within our scRNA-seq data, non-cell-autonomous dSPN-DEGs in Foxp1 samples were enriched for GO categories such as neuron projection (Table S5). Since projections onto the GPi were not altered in Foxp1 mice, dSPN-DEGs specific to Foxp1 samples are most likely responsible for the altered dSPN projection patterns found within Foxp1 animals. We therefore examined the overlap of dSPN-DEGs within Foxp1 (cell autonomous) and Foxp1 samples (non-cell autonomous) (Figure 5H). dSPN-DEGs specific to Foxp1 samples that are involved in neuron projection include Akap5, Asic2, Kirrel3, Cdh8, and Cntn4 (Figure 5H). Interestingly, Kirrel3, Cdh8, and Cntn4 are also ASD-risk genes (Figure 5H). These findings suggest deletion of Foxp1 within iSPNs alters the gene expression profiles within both iSPNs and dSPNs important for proper striatal projection patterning.

Distinct Behavioral Deficits with Cell-Type-Specific Deletion of Foxp1

We hypothesized that severe reduction of iSPNs and altered projection patterns with deletion of Foxp1 from iSPNs would result in altered motor behaviors. We therefore first tested behaviors classically characterized as being governed by striatal circuits, such as motor learning and activity levels. To test motor learning, we used the accelerating rotarod assay and found that Foxp1 and Foxp1 mice had significant deficits at remaining on the accelerating beam compared to control and Foxp1 mice (Figure 6A). This phenotype was not due to differences in grip strength (Figures S6A and S6B) or gait abnormalities (Figures S6C–S6F). Foxp1 and Foxp1 mice were also hyperactive in the open field behavioral paradigm compared to control mice (Figure 6B); however, no difference was observed in novel-cage activity between genotypes (Figure S6G). There was no difference in time spent in the periphery versus the center of the open field between genotypes (Figure 6C), suggesting no changes in anxiety-like behavior.

Figure 6.

Foxp1 Regulates Behaviors via Distinct Striatal Circuits

(A) Latency to fall was measured on the accelerating rotarod. Foxp1 and Foxp1 mice exhibit significant deficits. Data are represented as mean ± SEM, n = 11 Foxp1; n = 17 Foxp1; n = 18 Foxp1; n = 12 Foxp1; *p < 0.05, **p < 0.005, and ***p < 0.0001, two-way ANOVA with Sidak’s multiple comparisons test.

(B and C) Mice were tested within the open field paradigm with velocity (B) and percent time spent in the periphery versus center (C) plotted. Foxp1 and Foxp1 mice had significant increase in activity with no difference in percent time spent in the periphery and center. Data are represented as mean ± SEM; n = 22 Foxp1; n = 14 Foxp1; n = 17 Foxp1; n = 4 Foxp1; ***p <0.0001, one-way ANOVA with Sidak’s multiple comparisons test.

(D–F) Neonatal isolation vocalizations were measured at P4, P7, and P10. (D) The number of isolation calls were significantly reduced in Foxp1D1 mice. (E) Mean frequency (kHz) of the isolation calls was significantly altered in Foxp1 mice and at P4 within Foxp1 animals. (F) The call slope or “structure” of the call was significantly altered over postnatal development in Foxp1 pups and specifically at P10 within Foxp1 pups. Data are represented as mean ± SEM; n = 71 Foxp1; n = 47 Foxp1; n = 36 Foxp1; n = 11 Foxp1; *p < 0.05, **p < 0.005, and ***p < 0.0001, two-way ANOVA with Sidak’s multiple comparisons test.

(G) Representative images of nests.

(H) Foxp1 and Foxp1 mice produced nests with significantly lower quality scores compared to Foxp1 and Foxp1 mice. Data are represented as mean ± SEM; n = 7 Foxp1; n = 4 Foxp1; n = 5 Foxp1; n = 5 Foxp1; **p < 0.005, one-way ANOVA with Sidak’s multiple comparisons test.

(I and J) Associative fear memory was assessed using the fear-conditioning (FC) paradigm. All Foxp1 cKO mice displays deficits in cued FC (I) shown as the percent of time spent freezing. Only Foxp1 and Foxp1 mice displayed deficits in contextual FC (J). Data are represented as mean ± SEM n = 23 Foxp1; n = 22 Foxp1; n = 11 Foxp1; n = 15 Foxp1; *p < 0.05, **p < 0.005, and ***p < 0.0001, two-way ANOVA with Dunnett’s multiple comparisons test. See also Figure S6.

Since genetic variants in FOXP1 are strongly associated with ASD, we next examined ASD-relevant social communication behaviors. Using a maternal separation paradigm, we recorded pup USVs at three postnatal time points (P4, P7, and P10). We found that Foxp1 mice produced significantly fewer calls with altered call slope compared to control pups (Figures 6D and 6E). In addition, Foxp1 pups had significantly lower pitch at P4, while Foxp1 mice exhibited deficits in pitch across all developmental time points (Figure 6F). No significant USV changes were measured solely in Foxp1 pups, and no deficits in pup righting reflex were observed at P4 (Figure S6H). We also tested nest-building behavior, an important communal behavior in rodents (Deacon, 2006; Silverman et al., 2010), and found that Foxp1 and Foxp1 mice produced low-quality nests compared to control and Foxp2 nests (Figures 6G and 6H).

Because individuals with FOXP1 mutations are frequently comorbid for intellectual disability (Meerschaut et al., 2017; Siper et al., 2017), we next assessed whether learning and memory circuits were altered using the cued and contextual fear-conditioning (FC) paradigm (Figures 6I and 6J). All Foxp1 cKO mice had significantly reduced freezing behavior during cued-evoked fear memory recall (Figure 6I); however, only Foxp1 and Foxp1 mice showed significant deficits in context-evoked fear memory (Figure 6J). While hippocampal and amygdala circuits are classically associated with fear conditioning, striatal D1 receptors are also important for mediating proper contextual fear conditioning in mice (Ikegami et al., 2014). We also found that striosome-matrix architecture was more severely disrupted over postnatal development in Foxp1 and Foxp1 adult animals compared to control and Foxp1 mice (Figure S6I).

DISCUSSION

In this study, we use single-cell transcriptomics to examine the molecular mechanisms underlying striatal neuronal specification by sequencing thousands of striatal cells across control and cell-type-specific Foxp1 conditional mouse models. We show that Foxp1 influences striatal development through cell-type-specific molecular pathways and describe the molecular, functional, and behavioral consequences of Foxp1 deletion within distinct striatal circuits (Figure 7). We found that Foxp1 is specifically required for the development of iSPNs that localize to both matrix and striosome compartments. A subpopulation of iSPNs that remain with deletion of Foxp1 localized primarily to the matrix or striosome-matrix border. Future work will help resolve the mechanism underlying the loss of iSPNs with deletion of Foxp1. We also found that Foxp1 maintains striosome-matrix architecture through pathways in both iSPNs and dSPNs and found the cell-typespecific targets regulated by Foxp1 within each cell type. Further studying these Foxp1 targets could help elucidate how striatal neurochemical compartments are maintained. While we found that Cre turned on in D1- and D2-Cre strains within a 24-h window of each other, the precise embryonic timing of Foxp1 deletion in each cell type remains undetermined. Therefore, potential differences in the timing of Foxp1 loss in dSPNs versus iSPNs could play a role in the cell-type-specific effects observed in our data.

Figure 7.

Summary of Cellular, Structural, Functional, and Behavioral Findings within Cell-Type-Specific Foxp1 Conditional Knockout Mice

Foxp1 mice have an increase in eSPN subpopulations, reduced striosomal area, no gross SPN projection deficits, and distinct behavioral deficits relevant to social communication behavior and contextual fear conditioning. Foxp1 mice have a marked decrease in iSPN and increase in eSPN subpopulations; reduced striosomal area with few striosomal iSPNs; dSPN-GPi and iSPNGPe projection deficits; and deficits in motor learning, activity, and cued fear conditioning.

Deletion of Foxp1 in iSPNs led to fewer iSPN projections onto the GPe and fewer dSPN collaterals onto the GPi. Whether this non-cell-autonomous effect is caused from the loss of iSPNs or other mechanisms regulated by Foxp1 will need to be further investigated. dSPNs and iSPNs are known to form inhibitory axon collaterals onto neighboring SPNs and modulate their excitability (Tunstall et al., 2002, Taverna et al., 2008; Tecuapetla et al., 2009). iSPNs and dSPNs also cooperate together to intermix within the striosome and matrix compartments (Tinterri et al., 2018). We not only found that manipulation of iSPNs led to functional changes of dSPNs but also captured a molecular snapshot of this inter-SPN communication, including differentially expressed ASD-risk genes involved in neuron projection such as Cntn4, Cdh8, and Kirrel3. These findings help further our understanding of how both cell types cooperate to regulate proper striatal circuit formation.

FOXP1 is among a subset of genes repeatedly and significantly linked to ASD (Iossifov et al., 2014; Stessman et al., 2017; Satterstrom et al., 2020). Individuals with FOXP1 mutations are diagnosed with ASD or ASD features, and all reported cases are comorbid with intellectual disability, gross motor delays, and/or selective language impairments (Meerschaut et al., 2017; Siper et al., 2017). Mice with iSPN-deletion of Foxp1 caused significant motor disruptions. Concordant with our data, mice with ablated iSPNs or mice with Ppp1r1b deletion from iSPNs were also hyperactive in the open field (Bateup et al., 2010; Durieux et al., 2009). Adult mice with induced ablation of D2 receptors displayed severe motor learning impairments on the accelerating rotarod (Bello et al., 2017). These data indicate that loss of iSPNs with deletion of Foxp1 lead to significant motor-learning and activity deficits.

Pup USVs measure affective state and social behavior in mice (Boulanger-Bertolus et al., 2017; Silverman et al., 2010) and peak between P4 and P10 (Araujo et al., 2015). Disruption of neonatal call number and structure with deletion of Foxp1 within dSPNs is particularly interesting given the high coexpression of both Foxp1 and Foxp2 within this cell type and the ability of Foxp1 and Foxp2 to heterodimerize to regulate gene expression (Li et al., 2004). Foxp2 plays a critical role in the vocal behavior across many species, including humans, mice, and songbirds (Konopka and Roberts, 2016). We show that Cntnap2, a known shared target of Foxp1 and Foxp2 (O’Roak et al., 2011; Vernes et al., 2008), is significantly upregulated within dSPNs of both Foxp1 and Foxp1 samples. Variants in CNTNAP2 are also associated with ASD (O’Roak et al., 2011), and Cntnap2 KO mice have deficits in pup USV production and nest-building behavior (Peñ agarikano et al., 2011). We previously found that Foxp1 heterozygous mice display altered USV phenotypes, including deficits in call number, call structure, and pitch (Araujo et al., 2015). Since both D1- and D2-Cre strains express Cre in certain areas outside of the striatum (Figures S1H and S1I), we cannot completely rule out the possibility that loss of Foxp1 in other cell types contributes to our behavioral findings. However, recent studies characterizing mice with Foxp1 deletion in cortical-hippocampal circuits (Emx1; Foxp1 cKO mice) allows us to compare and contrast behaviors across mouse models with Foxp1 deleted in distinct regions/cell types (Araujo et al., 2017; Usui et al., 2017). Emx1; Foxp1 cKO mice produced fewer USVs, though no changes were observed in call structure or pitch (Usui et al., 2017). Here, we observed changes in all three parameters within Foxp1 and Foxp1 mice suggesting that Foxp1 regulates distinct aspects of mouse vocal behavior largely through cortical-striatal circuitry. Emx1; Foxp1 cKO mice also displayed deficits in nest-building behavior similar to deficits seen in Foxp1 and Foxp1 animals. These findings suggest that Foxp1 may regulate more general social-communication behaviors via cortico-striatonigral pathways.

Striosome compartments are smaller and architecturally disorganized with deletion of Foxp1 in iSPNs and/or dSPNs in the early postnatal striatum. Loss of striosome-matrix compartmentalization is particularly striking in adulthood with dSPN-specific deletion of Foxp1. These findings indicate that dSPN-targets regulated by Foxp1 exert a stronger influence over maintaining striatal neurochemical organization. Behaviors specific to Foxp1 mice include deficits in contextual fear memory recall, a known limbic-circuitry associated behavior. Striosomes receive preferential inputs from limbic subcortical regions, including the amygdala and bed nucleus of the stria terminalis (Smith et al., 2016); thus, inputs from these limbic regions targeting striosomes may be disrupted and contribute to the limbicassociated behavioral deficits seen in Foxp1 and Foxp1 mice. Additionally, Emx1; Foxp1 cKO mice did not show deficits in cued or contextual fear conditioning (Araujo et al., 2017). Therefore, Foxp1 is likely mediating fear-conditioned behaviors via specific disruption of striatal circuits.

While ASD is a genetically complex disorder, several studies have shown that striatal SPNs may be particularly vulnerable to ASD-linked mutations (Chang et al., 2015; Coe et al., 2019; Takata et al., 2018; Turner et al., 2017; Xu et al., 2014). Our study uncovers the molecular targets of Foxp1 in SPN subtypes and finds that Foxp1 regulates ASD-relevant behaviors via distinct striatal circuits. We show that iSPNs are particularly vulnerable with loss of Foxp1 and that Foxp1-regulated iSPN-targets are enriched for high-confidence ASD risk-genes, suggesting that striatopallidal circuitry might be particularly at risk with loss-offunction FOXP1 mutations. Our data provide important molecular insights for the development of future therapies targeting striatal circuits.

STAR★METHODS

LEAD CONTACT AND MATERIALS AVAILABILITY

All material and resource requests should be directed to and will be fulfilled by the Lead Contact, Genevieve Konopka (Genevieve. Konopka@utsouthwestern.edu). This study did not generate new unique reagents.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

All experiments were performed according to procedures approved by the UT Southwestern Institutional Animal Care and Use Committee (IACUCC# 2016–101-825). Foxp1 mice (Zhang et al., 2010) were provided by Dr. Haley Tucker and backcrossed to C57BL/6J for at least 10 generations to obtain congenic animals as previously described (Araujo et al., 2017; Usui et al., 2017). Drd1a-Cre (262Gsat, 030989-UCD) and Drd2-Cre (ER44Gsat, 032108-UCD) mice were obtained from MMRC. Drd2-eGFP (Gong et al., 2007) and Drd1-tdTomato (Ade et al., 2011) mice were provided by Dr. Craig Powell. We bred individual Cre or reporter lines to Foxp1 mice to obtain all Foxp1 cKO mice in one litter that were heterozygous for Cre or reporter transgene. Both male and female mice were used in all experiments. Mice used for single-cell RNA-sequencing and behavior experiments were not crossed with Drd1- or Drd2-reporter mice. Reporter mice were crossed with Foxp1 cKO lines for immunohistochemistry experiments and neuronal projection quantification. For timed-breeding experiments, we used the R26R-EYFP strain from Jackson Laboratory (stock no: 006148) bred to either Drd2-Cre; or Drd1-Cre; Foxp1 or Foxp1. Mice were maintained on a 12-hr light on/off schedule.

METHOD DETAILS

Protein isolation and immunoblotting

Striatal tissue was dissected, flash frozen, and stored at −80C before protein extraction. Protein was extracted from tissue using 1X RIPA Buffer (750mM NaCl, 250mM Tris-HCl pH7.4, 0.5% SDS, 5% Igepal, 2.5% Sodium deoxycholate, 5mM EDTA, 5mM NaVO4) with fresh protease inhibitor cocktail (10ul/ml), 10ul/ml of 100mM PMSF, and 25ul/ml of 200mM sodium orthovanadate. Tissue was homogenized in RIPA buffer using the TissueLyser LT (QIAGEN) with a sterile, stainless-steel bead for 1min at 50 Hz. Samples were agitated for 1hr at 4C, spun down at 12,000rpm for 15 min, and supernatant was transfer to a fresh tube. Protein was quantified using a standard Bradford assay (Bio-Rad) and 20ug of protein per sample were run on 10% SDS-Page gels. PVDF membranes (Bio-Rad, 162–0177) were incubated in blocking solution (1% Skim milk in TBS with 0.1% Tween-20) for 30 min at room temperature (RT) and probed with primary antibodies overnight at 4C. Membranes were washed with TBS-T (TBS with 0.1% Tween-20) and incubated with appropriate, species-specific fluorescent secondary antibodies in blocking solution for 1hr at RT, and washed in TBS-T. Images were collected using the Odyssey infrared imaging system (LI-COR Biosciences).

RNA isolation and quantitative real-time PCR

RNA from fresh or flash frozen tissue was harvested using miRNAeasy kit guidelines. RNA was converted to cDNA using recommended guidelines from SSIII Superscript Kit (Invitrogen) and qRT-PCR was performed using the CFX384 Real-Time System (Bio-Rad).

Immunohistochemistry

For P7 or P9 mice, rapid decapitation was performed. Brains were extracted and dropped into ice-cold 1X PBS for 1min before transfer into 4% PFA overnight. Brains were then transferred to 30% sucrose with 0.01% sodium azide for 48 hours. 35um coronal slices were made using a SM2000 R sliding microtome (Leica) and free-floating sections were stored in 1X PBS with 0.01% sodium azide. Slices were washed with 1X TBS and incubated for 30min in 3M glycine in 0.4% Triton-X, TBS. Slices were incubated in primary antibodies overnight at 4C, washed, and incubated in secondary antibodies for 1hr at room temperature. Slices were washed then mounted onto slides and allowed to dry overnight. Sections were incubated in DAPI solution (600nM in PBS) on the slide for 5 minutes and washed 3X with 1X PBS. Sections were allowed to dry before mounting coverslips using Prolong Diamond Antifade Mountant.

For embryonic time points, dams were anesthetized with CO2 and embryos were extracted rapidly and placed in a Petri dish with ice-cold 1X PBS. Embryos were transferred to clean Petri dishes with ice-cold 1X PBS twice before drop fixing in 4% PFA overnight. Embryos were transferred to 30% sucrose with 0.01% sodium azide for 48hrs before mounting embryos in OCT media and freezing at −20C. Sections were made at 20um on a Leica CM1950 Cryostat and directly mounted onto slides and allowed to dry before processing as described above.

Imaging and Analysis

Images were collected using a Zeiss Confocal laser scanning microscope (LSM880) and all image quantification was performed using Fiji image processing package. For iSPN quantification, 20X z stack images of dorsolateral, dorsomedial, and ventral striatum were taken within one hemisphere of four separate striatal sections from anterior to posterior per animal (3 images/section, 4 sections/animal, at least 3 animals/genotype). All images were taken within approximately similar sections across samples. D2-eGFP+ cells were manually quantified from maximum projection images within a 1024×1024 pixel field of view across all images and averaged per section. For striosome quantification, 10X z stack images were taken from one hemisphere of four separate striatal section from anterior to posterior per animal (4 sections/animal, at least 3 animals/genotype). Individual MOR+ patches were numbered, and area measurements summed for the total striosomal area measurement per section. Total striatal area was also measured per section to calculate the percentage of striosome area to total area per section. For iSPN within striosome compartment quantification, D2-eGFP+ cells located within MOR+ patches above the anterior commissure for each section (from images described above) were manually quantified. All D2R+ cells were automatically quantified within the same sections by using the ‘Analyze particles’ tool in Fiji with size set to 7.25-infinity (um2) (2–3sections/animal, 2–3animals/genotype). Differences between genotypes were assessed using a one-way ANOVA with multiple comparisons. For DEG quantification via IHC, 63X oil images were taken of each sample and target DEG corrected total cell fluorescence [CTCF = Integrated density − (Area of cell * average mean background measurements)] was measured within either D2R+ cells (eGFP+) or D1R+ (tdTom+) cells (15–40 cells/cell-type, 2–3 animals/genotype).

Antibodies

The following primary antibodies were used for either immunoblots (IB) or immunohistochemistry (IHC) experiments: chicken antiGFP (1:1,000, Aves Labs, GFP-1010), rabbit polyclonal anti-MOR (1:350, Millipore, AB5511), rabbit polyclonal anti-PDE1A (1:500, Proteintech, 12442–2-AP), rabbit polyclonal anti-PPP1R1B (1:1,000, Millipore, AB10518), goat anti-tdTomato (1:500, LifeSpan Biosciences, LS-C340696), mouse monoclonal anti-FOXP1 (1:500, Abcam, ab32010), rabbit polyclonal anti-FOXP1 (IHC:1:1,000, IB: 1:5,000 (Spiteri et al., 2007), rabbit polyclonal anti-Calbindin (1:500, Millipore AB1778 and 1:250, Swant CB-38a for DEG quantification), goat anti-FOXP2 (N-terminal) (1:500, Santa Cruz 21069), rabbit polyclonal anti-b-Tubulin (IB: 1:10,000, Abcam, ab243041), and mouse monoclonal anti-SOX4 (1:500, Abcam, ab243041). All IHC following secondary antibodies were used at a 1:1,000 dilutions Alexa Fluor 488 Donkey Anti-Chicken IgG (Thermo Fisher, 703–545-155), Alexa Fluor 555 Donkey Anti-Goat IgG (Thermo Fisher, A-21432), Alexa Fluor 647 Donkey Anti-Rabbit IgG (Thermo Fisher, 711–605-152), Alexa Fluor 647 Donkey Anti-Mouse IgG (Thermo Fisher, A-31571). For IB, the following secondary antibodies were used at a 1:10,000 dilution: IRDye 800CW Donkey anti-Rabbit IgG (Li-Cor, 925–32213) and IRDye 680RD Donkey anti-Rabbit IgG (Li-Cor, 925–68071).

Tissue processing for single-cell RNA-sequencing (scRNA-seq)

Mice (P9) were sacrificed by rapid decapitation and brains were quickly removed and placed in ACSF (126mM NaCl, 20mM NaHCO3, 20mM D-Glucose, 3mM KCl, 1.25mM NaH2PO4, 2mM of CaCl2 and MgCL2 freshly added) bubbled with 95%O2 and 5%CO2. Coronal slices at 500um were made using a VF-200 Compresstome in ACSF and transferred to a recovery chamber at room temperature in ACSF with 50uM AP5, 20uM DNQX, and 100nM TTX (ACSF+cb) (Tasic et al., 2016). Striatal punches including both dorsal and ventral striatum were taken from these slices and incubated in 1mg/ml of pronase in ACSF+cb for 5min. Punches were washed with ACSF+ 0.04% BSA twice and gently dissociated into single-cell suspension using polished Pasteur pipettes with 600um, 300um, and 150um opening diameters, sequentially. Cells were centrifuged and washed twice, filtered through Flowmi Tip 40uM strainers, and resuspended with ACSF+ 0.04% BSA. Cell viability was quantified using the trypan blue exclusion method and cell concentration was adjusted for targeted sequencing of 10,000 cells/sample using the 10X Genomics Single Cell 3′ Reagent Kits v2 protocol to prepare libraries (Zheng et al., 2017). A total of 16 mice (4 mice/genotype, 2 males and 2 females per genotype) were processed for single-cell sequencing. Libraries were sequenced using the McDermott Sequencing Core at UT Southwestern.

Pre-processing of Sequencing Data

Raw sequencing data was acquired from the McDermott Sequencing Core at UT Southwestern in the form of binary base call (BCL) files. BCL files were then de-multiplexed with the 10X Genomics i7 index (used during library preparation) using Illumina’s bcl2fastq v2.17.1.14 (Andrews, 2010) and mkfastq command from 10X Genomics CellRanger v2.1.1 tools (Zheng et al., 2017). Extracted paired-end fastq files (26 bp long R1 - cell barcode and UMI sequence information, 124 bp long R2 - transcript sequence information) were checked for read quality using FASTQC v0.11.5(Andrews, 2010). R1 reads were then used to estimate and identify real cells using whitelist command from UMI-tools v0.5.4 (Smith et al., 2017) program. A whitelist of cell-barcodes (putative real cells) and R2 fastq files were later used to extract reads corresponding to real cells only (excluding sequence information representing empty beads, doublets, low quality/degrading cells, etc.) using extract command from UMI-tools v0.5.4 (Smith et al., 2017). This step also appends the cell-barcode and UMI sequence information from R1 to read names in R2 fastq file. Extracted R2 reads were then aligned to reference mouse genome (MM10/GRCm38p6) from UCSC genome browser (Kent et al., 2002) and reference mouse annotation (Gencode vM17) using STAR aligner v2.5.2b (Dobin et al., 2013) allowing up to 5 mismatches. Uniquely mapped reads were then assigned to exons using featureCounts program from Subread package (v1.6.2) (Liao et al., 2014). Assigned reads sorted and indexed using Samtools v1.6 (Li et al., 2009) were then used to generate raw expression UMI count tables using count command from UMI-tools v0.5.4 (Smith et al., 2017) program. This raw expression matrix contains cells as rows and genes as columns and can be further used for downstream analysis such as normalization, clustering, differentially expressed genes, etc.

Clustering Analysis

Raw single-cell RNA-seq UMI count data was used for clustering analysis using Seurat R analysis pipeline (Butler et al., 2018). First, cells with more than 50,000 molecules (nUMI per cell) and cells with more than 10% mitochondrial content were filtered out to discard potential doublets and degrading cells. Also, genes from mitochondrial chromosome and chromosomes X and Y were removed as samples were from mixed genders. This dataset is referred to as primary filtered dataset. Post filtering, the raw UMI counts from primary filtered dataset were used for log-normalization and scaled using a factor of 10,000 and regressed to covariates such as number of UMI per cells and percent mitochondrial content per cell as described in Seurat analysis pipeline (Butler et al., 2018). To further identify the top variable genes, the data were used to calculate principal components (PCs). Using Jackstraw analysis, statistically significant PCs were used to identify clusters within the data using original Louvain algorithm as described in Seurat analysis pipeline followed by visualizing the clusters with uniform manifold approximation and projection (UMAP) in two dimensions (Becht et al., 2018). Genes enriched in each cluster compared to the remainder of the cells (adj. p value % 0.05 and log fold change > = 0.3) were identified as described in Seurat analysis pipeline. Genes corresponding to each cluster were used to identify the cell-type by correlating to genes expressed in previously published adult mouse striatal single cell data (Saunders et al., 2018). Cell-types were assigned to clusters based on (i) statistically significant enrichment of gene sets using the hypergeometric test (with a background of 7,500 genes, the number of expressed genes within our dataset) and (ii) expression weighted cell-type enrichment (EWCE) analysis (Skene and Grant, 2016) (https://github.com/NathanSkene/EWCE). Clusters that overlapped significantly with multiple cell-types were called for the most significant overlap (smallest Adj. P value) and analyzed for expression of top marker genes of known cell-types. Cells from clusters that fell into neuronal categories (referred to as secondary neuronal dataset) were used to re-cluster the cells to define specific spiny projection neuronal sub-types using a similar approach as described above. Note that two small clusters (Clusters-21, 22) that corresponded to excitatory cortical neurons and a cluster with less than 30 cells total (Cluster-24) were excluded from the secondary neuronal dataset UMAP plots to focus on striatal cell-types.

Differential Gene Expression (DEG) Analyses

Pairwise DEG analysis SPNs

For the spiny projection neuronal sub-type clusters identified using secondary neuronal dataset, pairwise differential gene expression analysis tests were performed within each cluster-pair using a Poisson likelihood ratio test from the Seurat R analysis pipeline (Butler et al., 2018) to identify genes enriched (adj. p value ≤ 0.05, |log2FC| > 0.25) in SPN sub-types.

Pseudobulk DEG analysis

Within the secondary neuronal dataset, neurons identified as either dSPNs (Drd1+) or iSPNs (Drd2+) regardless of cluster were combined into pools of cells segregated by genotypes. Differential expression within pools of dSPN or iSPNs of Foxp1 cKO samples were then compared to control samples using Poisson likelihood ratio test from the Seurat R analysis pipeline accounting for averaged expression differences in either dSPNs or iSPNs across genotypes irrespective of the identified clusters. Significant expression changes (adj. p value ≤ 0.05, |log2FC| > 0.3) reflected the differences in expression of genes in one specific cell population (dSPNs or iSPNs) across genotypes instead of detected clusters.

Down-sampled Dataset Analysis

Cells from the primary filtered dataset were used to randomly select the cells from each genotype matching the number of cells present in each genotype with the lowest representation of the cells (Foxp1 = 14466 cells, Foxp1 = 16,961 cells, Foxp1 = 9,898 cells, Foxp1 = 21,453 cells, using random sampling, the same number of cells from Foxp1, Foxp1and Foxp1 were matched to Foxp1). This is referred to as the primary down-sampled dataset. This dataset was further used to separate the cells into clusters and identify cell-types as described in the clustering analysis section above. Clusters corresponding to SPNs from the primary down-sampled dataset (referred to as the secondary down-sampled neuronal dataset) were re-clustered to identify SPN subtypes in a similar manner as described in the clustering section above.

TissueCyte Imaging and Quantification

STPT and image acquisition

Serial two-photon tomography (STPT) (Ragan et al., 2012), in which automated block face imaging of the brain is repetitively alternated with vibratome sectioning, was conducted on the TissueCyte 1000 platform using the manufacturer’s custom software for operation (Orchestrator). Mouse brains were perfusion-fixed in 4% paraformaldehyde and embedded in low-melting point oxidized agarose (4.5% w/v; Sigma #A0169). Vibratome sections were prepared at 75 μm thickness using a frequency of 70 Hz and a speed of 0.5 mm/sec. 185–190 total sections were collected of each brain. A 9 by 13 mosaic of tile images was collected at each level using lateral resolution of 0.875 μm/pixel. Optical sectioning was used to collect three z-planes within each 75 μm physical section to obtain 25 μm axial resolution. The two-photon excitation laser (Spectra Physics MaiTai DeepSee) was tuned to 920 nm to excite both eGFP and tdTomato. The emission fluorescence from the red, green and blue channels was independently collected using photomultiplier tube detectors. The tile images were saved to network attached servers and automatically processed to perform flat field correction and then stitched into single-channel 2D coronal sections in 16-bit .tif format using the manufacturer’s custom software (AutoStitcher).

Sample preparation and details

Mice (8–10 weeks, both male and female) were perfused with PBS followed by 4% PFA. Brains were removed and post-fixed overnight in 4% PFA at 4C. Samples were transferred to PBS + 0.1% sodium azide and stored at 4C until imaging. A total of 19 whole mouse brain images were collected in three cohorts for machine learning analysis according to their patterns of fluorophore expression. The first cohort consisted of 8 samples expressing tdTomato (detected predominantly in the red channel), the second cohort had 8 samples that expressed eGFP (detected predominantly in the green channel) and the third cohort consisted of 3 dual-labeled (eGFP + tdTomato) samples.

TissueCyte image processing and registration

STPT image processing was performed via BioHPC, an advanced computing cluster at UT Southwestern. All channels of the coronal sections were downsampled to 10 μm lateral resolution, intensity adjusted to fill the 16-bit range, and combined to form 3D image stacks using custom MATLAB software. The image stacks were then processed through a 3D median filter to remove high-contrast noise. The 3D image stacks were registered to Allen Institute for Brain Science Common Coordinate Framework (version 3, CCFv3) at 10 μm × 10 μm × 100 μm resolution using NiftyReg software (Modat et al., 2014). Briefly, registration involved three steps: (i) Affine transformation (reg-aladin) for global registration (ii) Cubic B-spline transformation (reg-f3d) to achieve local transformation and (iii) Resampling the transformed brains to Atlas coordinates (reg-resample). Registration transformations were established based on the red channel, then applied equally to all other data channels, including the probability maps (described below).

Interactive Image training for classifying signals of interest

The three raw channels of the 2D stitched coronal sections were downsampled to 1.5 μm lateral resolution. A maximum intensity projection of the three optical sections was produced for each physical section across all 3 color channels, creating an RGB image stack with the same number of 2D frames as physical sections (e.g., 185 or 190). Ilastik (Interactive learning and segmentation toolkit) (Sommer et al., 2011) software was deployed on BioHPC and used to train a pixel-wise random forest classifier to identify features of interest (e.g., fluorescent neuronal cell bodies and axonal projections). Three or four representative sections were chosen from the 185–190 image stack for model training. A supervised random forest model was trained by users to classify fluorescent features of interest (e.g., eGFP and/or tdTomato), and to distinguish them from other image features (e.g., bright microbubbles, empty space, autofluorescence) using the interactive features in Ilastik. An independent random forest model was trained for each of the image batches described above. The random forest classifiers were used to detect features of interest in all image sections, creating a “probability map” for each voxel in each 3D whole brain image. In these probability map images, the value of each voxel in each virtual channel (corresponding to each image feature, e.g., eGFP) represents the probability that the voxel includes information for the desired feature. These exported probability maps were registered to the CCFv3.0 using the transformation parameters using NiftyReg (reg-aladin).

Quantification and visualization

The features of interest in the registered probability maps were quantified by automatically segmenting brain regions of interest based upon CCFv3.0 volumetric annotations. Custom MATLAB software aggregated brain regions of interest (e.i., nucleus accumbens, caudate putamen, globus pallidus external and internal, substantia nigra pars compacta and pars reticulata), calculated the cumulative probabilities of all voxels in each region, and normalized these values by the volume of each structure. Cumulative probabilities normalized by the volume of each structure are normalized again to the highest signal intensity from a given region to obtain a value from 0–1 across all brain regions. This exported data matrix thus included normalized probability intensity values for each machine learning feature, each brain region of interest, and each brain. For visualization, the combined probability map stacks were rendered in 3D using the ClearVolume plugin for Fiji/ImageJ (Royer et al., 2015).

Behavior tests

Behavior analyses were performed on both male and female mice and no statistically significant differences were observed between sex. Subsequently, behavior data combines both male and female mice for downstream statistical analyses.

Open Field

Mice age 8–12 weeks were allowed to acclimate to the testing room for 1hr before being placed in a 55cm ×55cm × 36cm matrix (Phenome Technologies) and recorded for 30min. Total distance and velocity measurements were analyzed using Actimetrics LimeLight software.

Novel-cage activity

As previously described (Araujo et al., 2017), mice were moved into individual cages (18×28cm) with minimal bedding. Cage was placed into a dark Plexiglas box and the movements were measured using a Photobeam Activity System-Home Cage software for two hours. The number of beam breaks was recorded every 5 min and averaged over two hours for statistical analyses.

Rotarod

Following previously published methods (Araujo et al., 2015), mice (8–12 weeks) were acclimated to the testing room for 30min before placed in one lane of a 5-lane accelerating rotarod (Series 8 ITCC Life Science rotarod). The textured drum within the individual lanes was programed to accelerate from acceleration from 4–40 rpm within a maximum time frame of 300 s. Each mouse was positioned facing away from the experimenter. Latency to fall was recorded once the trial was initiated. Manual activation of the sensors occurred when an animal made a full rotation holding onto the drum. Animals received four trials per day (20min intervals) with lanes cleaned between animals with NPD over the course of three consecutive days.

Grip strength test

Grip strength was tested following previously published methods (Araujo et al., 2015). Briefly, following rotorad experiments, the forelimb and hindlimb grip strength mice were measured using Chatillon Force Measurement equipment. The forelimbs, followed by the hindlimbs, for each animal were tested first by placing forelimb paws on a mesh wire meter and pulling them away from the wire at constant force. Five consecutive measurements were recorded for both hindlimbs and forelimbs and averaged for a final grip strength measurement.

Nestlet behavior

Nesting behavior was analyzed using a previously published approach (Araujo et al., 2017; Deacon, 2006). Mice (8–12 weeks) were isolated into clean cages overnight with 3 g of intact nestlet. After 16–18 hr, the amount of unused nestlet was measured and images of the nests were taken to assess the quality and given a score.

Neonatal ultrasonic vocalization measurements

USVs were recorded as described previously (Araujo et al., 2015; 2017). Briefly, pups were isolated from dams at P4, P7, and P10 and placed into a soundproof container. USVs were recorded for 3min with an UltraSoundGate condenser microphone using Avisoft Bioacoustic software. Analysis of sound spectrograms was automatically performed using MATLAB codes (Rieger and Dougherty, 2016).

Digigait

Mice (8–12 weeks) were placed onto the transparent treadmill using the DigiGait Imaging System (Mouse Specifics, Inc) at 10 cm/sec. The speed was quickly increased to 20 cm/sec with a high-speed video camera mounted under the clear treadmill to capture images of all four paws at the 20 cm/sec speed. A section of video with at least 6–10 steps is analyzed and the paw placement is automatically detected and quantified by the software system. Right and left forelimb and hindlimb paw measurements were analyzed separately.

Fear Conditioning

Fear conditioning was measured using boxes with metal grid floors connected to a scrambled shock generator (Med Associates Inc., St. Albans). Mice were trained by placing them individually in the chamber for 2min before they received 3 tone-shock pairings (30sec white noise, 80dB tone, co-terminated with a 2 s, 0.5mA footshock, 1min intertrial interval). Twenty-four hours later, contextual memory was measured by placing the mice into the same chamber and measuring freezing behavior using the Med Associates software. Forty-eight hours post training, memory of the white noise cue was measured by placing mice in new environment, with altered floors, walls, different lighting, and a vanilla smell. Freezing was measured for 3 min and then noise cue was turned on for an additional 3 min and freezing was measured.

QUANTIFICATION AND STATISTICAL ANALYSIS

Statistical methods and code used for scRNA-seq and analysis are provided in the above methods sections. All statistical tests were performed using GraphPad Prism software (and p values obtained) for SPN projection analysis, behavior, and immunohistochemistry and are described in figure legends. No statistical methods were used to estimate sample size, but behavior cohorts were based on previously published papers (Araujo et al., 2015, 2017; Usui et al., 2017). Sample size for each experiment is indicated in figure legends.

DATA AND CODE AVAILABILITY

The sequencing data reported in this paper can be access at NCBI GEO with accession number GSE125290. Code that was used to perform data pre-processing, clustering and differential gene expression analysis is available at GitHub repository (https://github. “https://github.com/konopkalab/early-postnatal-striatal-single-cell-rna-seq”com/konopkalab/early-postnatal-striatal-single-cell-rna-seq “https://github.com/konopkalab/early-postnatal-striatal-single-cell-rna-seq“). All other acquired data are available upon request to G.K. (Genevieve.Konopka@utsouthwestern.edu).

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Chicken anti-GFP	Aves Labs	Cat#: GFP-1010 Lot#: GFP879484, RRID:AB_2307313
Rabbit polyclonal anti-MOR	Millipore	Cat#: AB5511 Lot#: 3131193, RRID:AB_177512
Rabbit polyclonal anti-PDE1A	Proteintech	Cat#:12442-2-AP, RRID:AB_2162893
Rabbit polyclonal anti-DARPP32	Millipore	Cat#: AB10518, RRID:AB_10807019
Mouse monoclonal anti-FOXP1 [JC12]	Abcam	Cat#: ab32010, RRID:AB_1141518
Rabbit polyclonal anti-FOXP1	(Spiteri et al., 2007)
Rabbit polyclonal anti-Calbindin	Millipore	Cat#: AB1778, RRID:AB_2068336
Rabbit polyclonal anti-Calbindin	Swant	Cat#: CB-38a, RRID:AB_2721225
Goat anti-FOXP2 (N-terminal)	Santa Cruz	Cat#: sc-21069, RRID:AB_2107124
Mouse anti-GAPDH	Millipore	Cat#: MAB374, RRID:AB_2107445
Mouse monoclonal anti-SOX4 [CL5634]	Abcam	Cat#: 243041
Alexa Fluor 488 Donkey Anti-Chicken IgG	Jackson ImmunoResearch Labs	Cat#: 703-545-155, RRID:AB_2340375
Alexa Fluor 555 Donkey Anti-Goat IgG	Thermo Fisher	Cat#: A-21432, RRID:AB_2535853
Alexa Fluor 647 Donkey Anti-Rabbit IgG	Jackson ImmunoResearch Labs	Cat#: 711-605-152, RRID:AB_2492288
Alexa Fluor 647 Donkey Anti-Mouse IgG	Thermo Fisher	Cat#: A-31571, RRID:AB_162542
IRDye 800CW Donkey anti-Rabbit IgG	Li-Cor	Cat#: 925-32213, RRID:AB_2715510
IRDye 680RD Donkey anti-Mouse IgG	Li-Cor	Cat#: 926-68072, RRID:AB_10953628
Chemicals, Peptides, and Recombinant Proteins
TTX	Tocris	Cat#: 1078
DNQX	Sigma-Aldrich	Cat#: D0540
APV	Tocris	Cat#: 3693
5X RIPA Buffer	In house
ProLong Diamond Antifade Mountant	Thermo Fisher	Cat#: P36970
DAPI	Thermo Fisher	Cat#: D1306
Dimethyl sulfoxide (DMSO)	Sigma	Cat#: D2438
Critical Commercial Assays
Chromium Single Cell 3′ Reagent Kits v2	10X Genomics, Inc	Cat#: 120237
miRNeasy mini kit	QIAGEN	Cat#:217004
SSIII First-strand super mix	Life Technologies	Cat#: 18080400
iTaq Universal SYBR Green Supermix	Bio-Rad	Cat#: 172-5124
Deposited Data
Raw and analyzed data	This paper	GEO: GSE125290
Adult striatal single-cell RNA-seq dataset	Saunders et al., 2018	GEO: GSE116470
Software and Algorithms
Bcl2fastq v2.17.1.14	Illumina Inc	https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html
CellRanger v2.1.1	10X Genomics	https://www.10xgenomics.com/solutions/single-cell/
FASTQC v0.11.5	Babraham Bioinformatics	https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
UMI Tools v0.5.4	(Smith et al., 2017)	https://github.com/CGATOxford/UMI-tools
STAR v2.5.2b	(Dobin et al., 2013)	https://github.com/alexdobin/STAR
Subread v1.6.2 (featureCounts)	WEHI	http://bioinf.wehi.edu.au/featureCounts/
Samtools v1.6	(Li et al., 2009)	https://github.com/samtools/samtools
Seurat v2.3.4	(Butler et al., 2018)	https://satijalab.org/seurat/, https://github.com/satijalab/seurat
EWCE	(Skene and Grant, 2016)	https://github.com/NathanSkene/EWCE
Toppgene	(Chen et al., 2009)	https://toppgene.cchmc.org
Ilastik version 1.2.2	(Sommer et al., 2011)	https://ilastik.org/
NiftyReg	(Modat et al., 2014)	https://github.com/KCL-BMEIS/niftyreg/wiki
Experimental Models: Organisms/Strains
Mouse: B6.FVB(Cg)-Tg(Drd1-cre)EY262Gsat	MMRRC	030989-UCD
Mouse: B6.FVB(Cg)-Tg(Drd2-cre)ER44Gsat	MMRRC	032108-UCD
Mouse: Drd1-tdTomato BAC Transgenic	(Ade et al., 2011)	N/A
Mouse: Drd2-eGFP BAC Transgenic	(Gong et al., 2007)	N/A
Oligonucleotides
Mus-Foxp1 F: CTACCGCTTCCATGGGAAAT	This paper	N/A
Mus-Foxp1 R: ACTGTGGTTGGCTGTTGTCA	This paper	N/A
Mus-Actin F: CCATCACAATGCCTGTGGTA	This paper	N/A
Mus-Actin R: CTAAGGCCAACCGTGAAAAG	This paper	N/A
Drd1-Cre Forward genotyping primers: GCTATGGAGATGCTCCTGATGGAA	MMRRC	N/A
Drd1-Cre Reverse genotyping primers: CGGCAAACGGACAGAAGCATT	MMRRC	N/A
Drd2-Cre Forward genotyping primers: GTGCGTCAGCATTTGGAGCAA	MMRRC	N/A
Drd2-Cre Reverse genotyping primers: CGGCAAACGGACAGAAGCATT	MMRRC	N/A
Drd1-tdTomato Forward genotyping primers: CTTCTGAGGCGGAAAGAACC	(Ade et al., 2011)	N/A
Drd1-tdTomato Reverse genotyping primers: TTTCTGATTGAGAGCATTCG	(Ade et al., 2011)	N/A
Drd2-eGFP Forward genotyping primers: CCCGAAGCTTCTCGAGGCGCGCCCTGTGCGTCAGCATTTGGAGCAAC	(Gong et al., 2007)	N/A
Drd2-eGFP Reverse genotyping primers: TCAGGGTCAGCTTGCCGTAGG	(Gong et al., 2007)	N/A
Foxp1-flox Forward genotyping primers: CCAGGGATCAGAGATTACTGTAGC	This paper	N/A
Foxp1-flox Reverse genotyping primers: CACCCTCTCCAAGTCTGCCTCAG	This paper	N/A