METHOD DETAILS

Nuclei were isolated from SC using detergent-based digestion and mechanical dissociation followed by a sucrose density gradient.¹⁹ Briefly, SC from E19, P4, P8, and P21 mice were manually isolated and homogenized in RNAase-free lysis buffer (0.32 M sucrose, 3 mM CaCl₂, 3 mM MgAc₂, 0.1 mM EDTA, 10 mM Tris-HCl,1 mM DTT, 0.1% Triton X-100 in DEPC-treated water) using a glass Dounce homogenizer on ice. For each sample, SC were collected from multiple mouse donors and pooled (i.e., two mice for E19 and five mice for all other time points). Care was taken to micro-dissect SC regions specifically without substantial contamination from deeper brain regions, particularly for the E19 brain in which the SC is less defined structurally. Animals of both sexes were used for all ages except for E19 which were all males. The homogenate was loaded into a polycarbonate ultracentrifuge tube containing sucrose solution (1.8 M sucrose, 3 mM MgAc₂, 1 mM DTT, 10 mM Tris-HCl in DEPC-treated water) in the bottom and centrifuged at 107,000 g for 2.5 h at 4°C. Supernatant was aspirated, and the nuclei containing pellet was incubated in RNAse-free 1x PBS, 0.04% BSA, 0.2 U/μL RNAse inhibitor on ice before resuspending the pellet. The nuclear suspension was filtered twice through a 30 μm cell strainer and counted using a Nexcelom Cellometer K2 before performing single-nucleus capture on the 10X Genomics 3′ v3 single cell RNA-Seq. Target capture of 2,000 nuclei per sample was used and the 10 × 3′ v3 scRNA-Seq library preparation performed and sequenced on the NovaSeq SP 100 (200,000 reads/nucleus) by the Oncogenomics Core Facility at University of Miami Miller School of Medicine.

After sequencing, Illumina output was processed using CellRanger v3.0.2. Base call files for each sample were demultiplexed. A pre-mRNA reference was generated with CellRanger mkref using the mm10 mouse genome. Each sample was aligned to the custom mm10 mouse reference genome using CellRanger. Sample reads were sequenced across two lanes and concatenated after alignment, resulting in a single count matrix per sample.

To distinguish nuclei-containing droplets from empty droplets, we performed cell calling on the unfiltered UMI count matrices using a combination of barcode-ranking and the empty-droplet detection emptyDrops function as implemented in the DropletUtils R package.⁵⁵ First, nuclei were ranked according to total UMI count and visualized in a log-total UMI vs. log rank plot (Figure S1A). A spline curve was fit to the data to identify “knee” and inflection points, and cells with total UMI count above the knee were considered nuclei-containing droplets. Next, we used the emptyDrops algorithm to further distinguish empty droplets from nuclei for nuclei with lower total UMI counts.

To remove potential doublets, we applied the Python package Scrublet⁶⁰ to each individual sample using default parameters. In brief, Scrublet simulates multiplets by sampling from the data and builds a nearest-neighbor-based classifier. Cells with high doublet scores were flagged and removed.

We performed further quality control based on metrics such as total UMI, number of unique genes detected, and mitochondrial transcript content (Figure S1B). Lower-bound thresholds for total UMI and unique gene detection rates were determined by computing three absolute median deviations (MADs) below the median. Across all samples, mean total UMI was 10,556 (E19, 7,457; P4, 9,882; P8, 10,957; P21, 13,928), mean total genes detected was 3,550 (E19, 3,030; P4, 3,574; P8, 3,643; P21, 3,953), and mean mitochondrial transcript content was 0.62% (E19, 0.70%; P4, 0.98%; P8, 0.41%; P21, 0.39%). Remaining high-quality nuclei were used for downstream analysis.

Since donor cells were not tagged prior to pooling, we could not measure the relative cell yield per donor mouse. To approximate the distribution of sexes to cell-types and subtypes, we calculated the percentage of Xist+ cells per sample and per cluster. We found that samples E19, P4, P8, and P21 were 0.9%, 46.0%, 40.7%, and 68.1% Xist+ cells, respectively, with a mean of 37.7% across all cells from all samples. From this, we speculated that sex-specific gene expression may contribute to batch effects and/or sex-specific clusters. However, we also found that neuron subtypes ranged in Xist+ percentage from 15.6% to 51.6% with a mean of 37.5%. This suggests that sex-specific gene expression did not strongly drive neuron subtype clustering results and that our neuron subtype analysis captures gene expression variation irrespective of sex.

We first performed standard single-nucleus RNAseq analysis using the Seurat R package (v4.2.1).⁶¹ We observed significant batch effects between samples from each developmental time point; for example, neurons from P21 clustered separately from neurons from all other time points due to detection of contaminant ambient RNAs such as Plp1 (Figure S1C). Based on this observation, we determined that Data Integration was necessary.

To better identify shared and unique cell types across all developmental time points, we performed integrated analysis as outlined in Seurat’s Data Integration workflow.⁶² In brief, after UMI count matrices were log-normalized, the top 2,500 variable genes and first 15 principal components were used for dimensional reduction and clustering. Cell types were identified using a combination of DE testing, comparisons against reference data sets, and prior knowledge of cell type-specific marker genes. For DE testing to identify markers, we used the FindAllMarkers function in Seurat using default parameters. For comparisons against reference data sets, we used the SingleR R package.³⁵ See methods section on “Comparison of neuronal subtypes to reference data” for description of the SingleR method. After cell type identification, DEGs were recomputed again using the FindAllMarkers function. We principally used the following genes for cell type identification: Excitatory neurons, Slc17a6; Inhibitory neurons, GAD1; Astrocytes, Aqp4 and Gfap; Oligodendrocyte-lineage cells, Cspg4, Bmp4, Mbp; Dividing cells, Mki67; Microglia, P2ry12; Endothelial cells, Cldn5; Vascular Leptomeningeal cells, Col1a1; Epithelial cells, Cdh1.

To investigate neuronal heterogeneity and identify neuronal sub-types, we performed a similar integrated analysis as described above with certain modifications. First, we used the top 3000 variable genes and top 10 principal components for dimensional reduction and clustering. Principal components were determined using the “elbow” plot heuristic. We also set the “resolution” parameter in the FindClusters function to 0.72 based on empirical observations that at this resolution neuron clusters could be identified using single or a small set of genes. Neuronal subtypes were annotated using these parameters. Neuronal subtype marker genes were computed using the FindAllMarkers function and filtered by taking the top 2 genes by p value and then by log(fold-change). To better identify sub-type markers within each of the excitatory and inhibitory neuron classes, we extracted each class and reperformed DE tests via FindAllMarkers. We used the default Seurat log2(fold-change) threshold of 0.25 and adjusted p value threshold of 0.05 for determining the number of DEGs (Figures S2E and S2F). We further quantified the similarity between neuron subtypes by constructing a dendrogram relating the average expression profile of neuron subtypes using the same genes used for cluster analysis (Figure S2D).

We also sought to identify gene expression changes in all excitatory or inhibitory neurons across development using the FindMarkers function. We observed that many of the top DEGs between time points were non-neuronal, e.g., the gene Ttr, which has been shown to be specific to cells of the choroid plexus in the CNS.⁶³ We reasoned that many of these genes may be derived from the ambient RNA during droplet processing.⁶⁴ To mitigate the effects of ambient RNA contamination in DE testing, we applied ambient profile estimation algorithms in the DropletUtils R package (Figure S1A). In brief, unfiltered UMI count matrices, containing transcript quantifications for all 10X Chromium lipid encapsulations, were used to estimate the ambient RNA profile based on expression data from low total UMI droplets. We used the ambientContribMaximum function to then filter out genes from DE test results which had an average maximum ambient RNA count contribution greater than 20% across all samples in the comparison. This approach was also applied in identifying gene expression changes across developmental time points for non-neuronal populations.

To further characterize the global changes in neuronal gene expression across development, we combined all neuronal cells into a single group and performed time point comparisons as well as ambient RNA filtering (Figure S3A). We performed Gene Ontology enrichment analysis for Biological Process terms using these DEGs (Figure S3B) using the topGO R package version 2.50.0.

To investigate potential development trajectories between neuronal subtypes, we used the Monocle3 R package.⁵⁷ We ported the Seurat dataset into Monocle3 and followed the Monocle3 pipeline except for the following parameters: for preprocess_cds() we “num_dim” to 20 for the PCA; we performed batch correction using align_cds(); we set the root node using the graphical interface prompted by order_cells(). To identify DEGs along trajectories, we used the graph_test() function using default parameters. We clustered DEGs using the find_gene_modules() function and setting the resolution parameter to range from 10^−6 to 10^−1 with increasing exponent.

To compare the neuronal subtypes identified in our study to neuronal subtypes described previously in other reports studying SC neurons, we applied the SingleR algorithm from the SingleR R package. In brief, SingleR first identifies marker genes for each neuron subtype label in the reference data in a pairwise manner. These genes are then used to compute Spearman correlations between the gene expression profiles of cells from the query dataset and cells of the neuron subtypes from the reference dataset. For each query cell, a per-subtype distribution is generated from the correlation values against cells from that neuron subtype label. For that query cell, the per-subtype label is defined as a fixed quantile (default 0.8) of this distribution. The subtype label with the highest score becomes the predicted subtype of the query cell. Heatmaps demonstrating the per-reference-label contribution to each neuronal subtype in the current study were generated using these labels.

To perform an integrated analysis, we pulled data from Cheung et al., Tsai et al., and the current study. For each dataset, we first performed an initial cluster analysis to identify neuronal cells using the marker genes Slc17a6, GAD1, and GAD2 and took the subset of data along genes which were shared among all datasets. We then computed the top variable genes for dimensional reduction using the modeleGeneVar() function from the BioConductor package scran.⁵⁸ We set the “block” argument to individually sequenced samples as our initial cluster analyses revealed significant batch effects between and within studies. We performed batch correction using the Seurat integration pipeline using default parameters and set the PCA dimensions to 20 for dimensional reductions.

To perform a comprehensive query of cell adhesion and axon guidance molecules, we pulled gene sets from the following sources: extracellular matrix and adhesion molecules were pulled from GeneCopoeia’s ExProfile Extracellular Matrix and Adhesion Molecules gene panel; axon guidance molecules were pulled from the KEGG pathway database using the pathway ID “mmu04360”. To investigate the heterogeneity of expression of these molecules in SC nuclei, we performed integrated cluster analysis using the Seurat integration pipeline as described above using these gene sets.

For Cre-dependent anterograde labeling of Pax7 expressing neurons in the superficial SC (sSC), approximately 100 nL of AAV2–CAG–FLEX–GFP (University of North Carolina Vector Core) was injected into the sSC of Pax7-Cre mice (10 weeks old). Injection coordinates were as follows: posterior from bregma, lateral from midline, and depth in mm, 4.16, 0.2, 0.5, and 1.0–1.2, respectively. Three to 4 weeks after AAV injection, mice were anesthetized and then transcardially perfused with 4% paraformaldehyde in phosphate buffered saline (PBS). Brains were dissected and postfixed in 4% paraformaldehyde in PBS for 16 h, and cryoprotected in 30% sucrose in PBS for 2–3 days. Brains were embedded in OCT compound (Tissue-Tek) and coronal sections (20 μm) were cut using a cryostat. Sections were immunostained by incubating in primary antibodies in 5% Normal Goat Serum in PBS with 0.3% Triton X-overnight at 4°C. Primary antibodies used were: RFP (Rockland 600–401-379S, 1:1000) and GFP (Abcam ab13970, 1:2000). Following primary antibody incubation, sections were washed and incubated in species-appropriate Alexa Fluor IgG (H + L) secondary antibodies (Invitrogen, 1:500) at room temperature for 1 h. Slides were mounted using Vectashield with DAPI (Vector Laboratories H-1200). Images were obtained using a Nikon Eclipse Ti fluorescent microscope or an Olympus FluoView 1000 confocal microscope.

RNAscope FISH was performed on 20 μm thickness coronal brain sections from adult mice (8 weeks old) using the RNAscope Multiplex Fluorescent v2 Assay (ACD Biotechne, Catalog No. 323100) according to the manufacture’s protocol. Target probes used are listed in the key resources table. TSA-based fluorophores were from PerkinElmer (TSA Plus Fluorescein, PN NEL741001KT; TSA Plus Cyanine 3, PN NEL744001KT; TSA Plus Cyanine 5, NEL745001KT). Images were acquired using an Olympus Confocal FV1000 microscope or an Andor Dragonfly confocal microscope.