A transcriptomic atlas of mouse cerebellar cortex comprehensively defines cell types

Animals

Nuclei suspensions for mouse (C57BL/6J, Jackson Labs) cerebellum profiles were generated from 2 female and 4 male adult mice (60 days old), 1 male E18 mouse, 1 male P0 (newborn) mouse, 1 female P4 (4 days old) mouse, 1 female P8, 2 male P12 and 2 female P16 mice. Adult mice were group-housed with a 12-h light-dark schedule and allowed to acclimate to their housing environment for two weeks after arrival. Timed pregnant mice were received and euthanized to yield E18 mice 6 days after arrival. Newborn mice were housed as individual litters for up to 16 days. All experiments were approved by and in accordance with Broad IACUC protocol number 012-09-16.

Brain preparation

At E18, P0, P4, P8, P12, P16 and P60, C57BL/6J mice were anaesthetized by administration of isoflurane in a gas chamber flowing 3% isoflurane for 1 min. Anaesthesia was confirmed by checking for a negative tail and paw pinch response. Mice were moved to a dissection tray and anaesthesia was prolonged via a nose cone flowing 3% isoflurane for the duration of the procedure. Transcardial perfusions were performed on adult, pregnant (E18), P8, P12 and P16 mice with ice-cold pH 7.4 HEPES buffer containing 110 mM NaCl, 10 mM HEPES, 25 mM glucose, 75 mM sucrose, 7.5 mM MgCl₂, and 2.5 mM KCl to remove blood from the brain. P0 and P4 mice were unperfused. The brain was removed from P60, P8, P12 and P16 mice and frozen for 3 min in liquid nitrogen vapour. E18, P0 and P4 mice were sagittally bisected after similarly freezing their brains in situ. All tissue was moved to −80 °C for long-term storage. A detailed protocol is available at protocols.io (https://doi.org/10.17504/protocols.io.bcbrism6).

Generation of cerebellar nuclei profiles

Frozen adult mouse brains were securely mounted by the frontal cortex onto cryostat chucks with OCT embedding compound such that the entire posterior half including the cerebellum and brainstem were left exposed and thermally unperturbed. Dissection of each of 16 cerebellar vermal and cortical lobules was performed by hand in the cryostat using an ophthalmic microscalpel (Feather safety Razor P-715) pre-cooled to −20 °C and donning four surgical loupes. Whole E18, P0, P4, P8, P12 and P16 mouse cerebella were similarly curated by dissecting rhombomeric cerebellar rudiments from sagittal frozen brain hemispheres using a pre-cooled 1-mm disposable biopsy punch (Integra Miltex). Each excised tissue dissectate was placed into a pre-cooled 0.25 ml PCR tube using pre-cooled forceps and stored at −80 °C. Nuclei were extracted from this frozen tissue using gentle, detergent-based dissociation, according to a protocol available at protocols.io (https://doi.org/10.17504/protocols.io.bck6iuze) adapted from one provided by the McCarroll laboratory (Harvard Medical School), and loaded into the 10x Chromium V3 system. Reverse transcription and library generation were performed according to the manufacturer’s protocol.

Floating slice hybridization chain reaction on acute slices

Acute cerebellar slices containing Alexa 594-filled patched cells were fixed as described and stored in 70% ethanol at 4 °C until hybridization chain reaction (HCR). They were then subjected to a ‘floating slice HCR’ protocol in which the recorded cells could be simultaneously re-imaged in conjunction with HCR expression analysis in situ and catalogued as to their positions in the cerebellum. A detailed protocol (https://doi.org/10.17504/protocols.io.bck7iuzn) was performed using the following HCR probes and matching hairpins purchased from Molecular Instruments: glutamate metabotropic receptor 8 (Grm8) lot number PRC005, connexin 36 (Gjd2) lot number PRD854 and PRA673, cadherin22 (Cdh22) lot number PRC011, neurexophilin 1 (Nxph1) lot number PRC675 and PRC466, leucine-rich glioma-inactivated protein 2 (Lgi2) lot number PRC012, somatostatin (Sst) lot number PRA213 and sortilin related VPS10 domain containing receptor 3 (Sorcs3) lot number PRC004. Amplification hairpins used were type B1, B2 and B3 in 488 nm, 647 nm and 546 nm respectively.

Patch fill and HCR co-imaging

After floating-slice HCR, slices were mounted between no.1 coverslips with antifade compound (ProLong Glass, Invitrogen) and images were collected on an Andor CSU-X spinning disk confocal system coupled to a Nikon Eclipse Ti microscope equipped with an Andor iKon-M camera. The images were acquired with an oil immersion objective at 60×. The Alexa 594 patched cell backfill channel (561 nm) plus associated HCR probe/hairpin channels (488 nm and 647 nm) were projected through a 10–20-μm thick z-series so that an unambiguous determination of the association between the patch-filled cell and its HCR gene expression could be made. Images were processed using Nikon NIS Elements 4.4 and Nikon NIS AR.

Human brain and nuclei processing

Human donor tissue was supplied by the Human Brain and Spinal Fluid Resource Center at UCLA, through the NIH NeuroBioBank. This work was determined by the Office of Research Subjects Protection at the Broad Institute not to meet the definition of human subjects research (project ID NHSR-4235).

Nuclei suspensions from human cerebellum were generated from two neuropathologically normal control cases—one female tissue donor, aged 35, and one male tissue donor, aged 36. These fresh frozen tissues had post mortem intervals of 12 and 13.5 h respectively, and were provided as whole cerebella cut into four coronal slabs. A sub-dissection of frozen cerebellar lobules was performed on dry ice just before 10x processing and nuclei were extracted from this frozen tissue using gentle, detergent-based dissociation, according to a protocol available at protocols.io (https://doi.org/10.17504/protocols.io.bck6iuze).

Electrophysiology experiments

Acute parasagittal slices were prepared at 240-μm thickness from wild-type mice aged P30–P50. Mice were anaesthetized with an intraperitoneal injection of ketamine (10 mg kg⁻¹), perfused transcardially with an ice-cold solution containing (in mM): 110 choline chloride, 7 MgCl₂, 2.5 KCl, 1.25 NaH₂PO₄, 0.5 CaCl₂, 25 glucose, 11.5 sodium ascorbate, 3 sodium pyruvate, 25 NaHCO₃, 0.003 (R)-CPP, equilibrated with 95% O₂ and 5% CO₂. Slices were cut in the same solution and were then transferred to artificial cerebrospinal fluid (ACSF) containing (in mM) 125 NaCl, 26 NaHCO₃, 1.25 NaH₂PO₄, 2.5 KCl, 1 MgCl₂, 1.5 CaCl₂ and 25 glucose equilibrated with 95% O₂ and 5% CO₂ at approximately 34 °C for 30 min. Slices were then kept at room temperature until recording.

All UBC recordings were done at 34–36 °C with (in μM) 2 (R)-CPP, 5 NBQX, 1 strychnine, 10 SR95531 (gabazine) and 1.5 CGP in the bath to isolate metabotropic currents. Loose cell-attached recordings were made with ACSF-filled patch pipettes of 3–5 MΩ resistance. Whole-cell voltage-clamp recordings were performed while holding the cell at −70 mV with an internal containing (in mM): 140 KCl, 4 NaCl, 0.5 CaCl₂, 10 HEPES, 4 MgATP, 0.3 NaGTP, 5 EGTA 5, and 2 QX-314, pH adjusted to 7.2 with KOH. Brief puffs of glutamate (1 mM for 50 ms at 5 psi) were delivered using a Picospritzer II (General Valve Corp.) in both cell-attached and whole-cell configuration to assure consistent responses. The heat map of current traces from all cells are sorted by the score over the first principal axis after singular value decomposition (SVD) of recordings over all cells.

For whole-cell recordings with pharmacology, we used an K-methanesulfonate internal containing (in mM): 122 K-methanesulfonate, 9 NaCl, 9 HEPES, 0.036 CaCl₂, 1.62 MgCl₂, 4 MgATP, 0.3 GTP (Tris salt), 14 creatine phosphate (Tris salt), and 0.18 EGTA, pH 7.4. A junction potential of −8 mV was compensated for during recording. 300nM TTX was added to the ACSF in conjunction with the synaptic blockers listed above. Three pipettes filled with ACSF containing 1 mM glutamate, 100 μM DHPG or 1 μM LY354740 were positioned within 20 μm of the recorded cell. Pressure applications of each agonist were delivered at 10 psi with durations of 40–50 ms. Agonist applications were separated by 30 s. Two to three trials were collected for each agonist.

MLI recordings were performed at approximately 32 °C with an internal solution containing (in mM) 150 K-gluconate, 3 KCl, 10 HEPES, 3 MgATP, 0.5 GTP, 5 phosphocreatine-tris₂ and 5 phosphocreatine-Na₂, 2 mg ml⁻¹ biocytin and 0.1 Alexa 594 (pH adjusted to 7.2 with KOH, osmolality adjusted to 310 mOsm kg⁻¹). Visually guided whole-cell recordings were obtained with patch pipettes of around 4 MΩ resistance pulled from borosilicate capillary glass (BF150-86-10, Sutter Instrument). Electrophysiology data was acquired using a Multiclamp 700B amplifier (Axon Instruments), digitized at 20 kHz and filtered at 4 kHz. For isolating spikelets in MLI recordings, cells were held at −65 mV in voltage clamp and the following receptor antagonists were added to the solution (in μM) to block synaptic currents: 2 (R)-CPP, 5 NBQX, 1 strychnine, 10 SR95531 (gabazine), 1.5 CGP. All drugs were purchased from Abcam and Tocris. To obtain an input-output curve, MLIs were maintained at 60–65 mV with a constant hyperpolarizing current, and 250 ms current steps ranging from −30 pA to +100 pA were injected in 10 pA increments. To activate the hyperpolarization-evoked current (I_h), MLIs were held at −65 mV and a 30 pA hyperpolarizing current step of 500 ms duration was injected. The amplitude of I_h was calculated as the difference between the maximal current evoked by the hyperpolarizing current step and the average steady-state current at the end (480–500 ms) of the current step. Capacitance and input resistance (R_i) were determined using a 10 pA, 50 ms hyperpolarizing current step. To prevent excessive dialysis and to ensure successful detection of mRNAs in the recorded cells, the total duration of recordings did not exceed 10 min. Acquisition and analysis of electrophysiological data were performed using custom routines written in MATLAB (Mathworks), IgorPro (Wavemetrics), or AxoGraphX. Data are reported as median ± interquartile range, and statistical analysis was carried out using the Mann–Whitney or Fisher’s exact test, as indicated. Statistical significance was assumed at P < 0.05.

To determine the presence of spikelets, peak detection was used to generate event-triggered average waveforms with thresholds based on the mean absolute deviation (MAD) of the raw trace. Spikelet recordings were scored for the presence of spikelets blind to the molecular identity of the cells. The analysis was restricted to cells recorded in the presence of synaptic blockers.

Imaging and analysis

MLIs were filled with 100 μM Alexa-594 via patch pipette to visualize their morphology using two-photon imaging. After completion of the electrophysiological recordings the patch electrode was retracted slowly and the cell resealed. We used a custom-built two-photon laser-scanning microscope with a 40×, 0.8 numerical aperture (NA) objective (Olympus Optical) and a pulsed two-photon laser (Chameleon or MIRA 900, Coherent, 800 nm excitation). DIC images were acquired at the end of each experiment and locations of each cell within the slice were recorded. Two-photon images were further processed in ImageJ.

Tissue fixation of acute slices

After recording and imaging, cerebellar slices were transferred to a well-plate and submerged in 2–4% PFA in PBS (pH 7.4) and incubated overnight at 4 °C. Slices were then washed in PBS (3 × 5 min) and then kept in 70% ethanol in RNase-free water until HCR was performed.

Preprocessing of sequencing reads

Sequencing reads from mouse cerebellum experiments were demultiplexed and aligned to a mouse (mm10) premrna reference using CellRanger v3.0.2 with default settings. Digital gene expression matrices were generated with the CellRanger count function. Sequencing reads from human cerebellum experiments were demultiplexed and aligned to a human (hg19) premrna reference using the Drop-seq alignment workflow², which was also used to generate the downstream digital gene expression matrices.

Estimation of adequate rare cell type detection

To estimate the probability of sufficiently sampling rare cell types in the cerebellum as a function of total number of nuclei sampled, we used the approach proposed by the Satija laboratory (https://satijalab.org/howmanycells), with the assumption of at most 10 very rare cell types, each with a prevalence of 0.15%. We derived this minimum based on the observed prevalences of the two rarest cell types we identified (OPC_4, Purkinje_Aldoc_2). We set 70 cells as the threshold for sufficient sampling, and calculated the overall probability as a negative binomial (NB) density:

in which k = 70, P = 0.0015, m = 10, and n represents the total number of cells sampled.

Cell type clustering and annotation

After generation of digital gene expression matrices as described above, we filtered out nuclei with fewer than 500 UMIs. We then performed cell type annotation iteratively through a number of rounds of dimensionality reduction, clustering, and removal of putative doublets and cells with high mitochondrial expression. For the preliminary clustering step, we performed standard preprocessing (UMI normalization, highly variable gene selection, scaling) with Seurat v2.3.4 as previously described³⁸. We used principal component analysis (PCA) with 30 components and Louvain community detection with resolution 0.1 to identify major clusters (resulting in 34 clusters). At this stage, we merged several clusters (primarily granule cell clusters) based on shared expression of canonical cell type markers, and removed one cluster whose top differentially expressed genes were mitochondrial (resulting in 11 clusters).

For subsequent rounds of cluster annotation within these major cell type clusters, we applied a variation of the LIGER workflow previously described¹², using integrative non-negative matrix factorization (iNMF) to limit the effects of sample- and sex-specific gene expression. In brief, we normalized each cell by the number of UMIs, selected highly variable genes⁷ and spatially variable genes (see section below), performed iNMF, and clustered using Louvain community detection (omitting the quantile normalization step). Clusters whose top differentially expressed genes indicated contamination from a different cell type or high expression of mitochondrial genes were removed during the annotation process, and not included in subsequent rounds of annotation. This iterative annotation process was repeated until no contaminating clusters were identified in a round of clustering. Differential expression analysis within rounds of annotation was performed with the Wilcoxon rank sum test using Seurat’s FindAllMarkers function. Comprehensive differential expression analysis across all 46 final annotated clusters was performed using the wilcoxauc function from the presto package³⁹. A full set of parameters used in the LIGER annotation steps and further details can be found in Supplementary Table 4.

For visualization as in Fig. 1b, we merged all annotated high-quality nuclei and repeated preliminary preprocessing steps before performing UMAP using 25 principal components.

Integrated analysis of human and mouse data

After generation of digital gene expression matrices for the human nuclei profiles, we filtered out nuclei with fewer than 500 UMIs. We then performed a preliminary round of cell type annotation using the standard LIGER workflow (integrating across batches) to identify the primary human interneuron populations (UBCs, MLIs and PLIs, Golgi cells, granule cells; based on the same markers as in Supplementary Table 1). We repeated an iteration of the same workflow for the four cell populations specified above (with an additional quantile normalization step) in order to identify and remove putative doublet and artefactual populations. Finally, we performed iNMF metagene projection as previously described⁴⁰ to project the human datasets into latent spaces derived from the corresponding mouse cell type datasets. We then performed quantile normalization and Louvain clustering, assigning joint clusters based on the previously annotated mouse data clusters. For the granule cell joint analysis, we first limited the mouse data to include only the five cerebellar regions sampled in human data collection (lobules II, VII, VIII, IX and X). For the Golgi cell joint analysis, we performed iNMF (integrating across species), instead of metagene projection.

Spatially variable gene selection

To identify genes with high regional variance, we first computed the log of the index of dispersion (log variance-to-mean ratio, logVMR) for each gene, across each of the 16 lobular regions. Next, we simulated a Gaussian null distribution whose centre was the logVMR mode, found by performing a kernel density estimation of the logVMRs (using the density function in R, followed by the turnpoints function). The standard deviation of the Gaussian was computed by reflecting the values less than the mode across the centre. Genes whose logVMRs were in the upper tail with P < 0.01 (Benjamini–Hochberg adjusted) were ruled as spatially variable. For the granule cell and PC cluster analyses, adjusted P-value thresholds were set to 0.001 and 0.002, respectively.

Cluster regional composition testing and lobule enrichment

To determine whether the lobule composition of a cluster differs significantly from the corresponding outer level cell type lobule distribution, we used a multinomial test approximated by Pearson’s chi-squared test with k − 1 degrees of freedom, in which k was the total number of lobules sampled (16). The expected number of nuclei for a cluster i and lobule j was estimated as follows:

where N_i is the total number of nuclei in cluster i and N_j is the number of nuclei in lobule j (across all clusters in the outer level cell type, as defined below). The resulting P values were FDR-adjusted (Benjamini–Hochberg) using the p.adjust function in R.

Lobule enrichment (LE) scores for each cluster i and each lobule j were calculated by:

in which n_ij is the observed number of nuclei in cluster i and lobule j, and N_j is the number of nuclei in lobule j (across all clusters in the outer level cell type). For this analysis, we used coarse cell type definitions shown coloured in the Fig. 2a, and merged the PLI clusters. For lobule composition testing and replicate consistency analysis below, we downsampled granule cells to 60,000 nuclei (the next most numerous cell type were the MLI and PLI clusters with 45,767 nuclei).

To determine the consistency of lobule enrichment scores across replicates in each region, we designated two sets of replicates by assigning nuclei from the most represented replicate in each region and cluster analysis to ‘replicate 1’ and nuclei from the second most represented replicate in each region to ‘replicate 2’. This assignment was used because not all regions had representation from all individuals profiled, and some had representation from only two individuals. We calculated lobule enrichment scores for each cluster using each of the replicate sets separately; we then calculated the Pearson correlation between the two sets of lobule enrichment scores for each cluster. We would expect correlation to be high for clusters when lobule enrichment is biologically consistent. We note that one cluster (Purkinje_Aldoc_2), was excluded from the replicate consistency analysis as under this design, it had representation from only a single aggregated replicate. However, we confirmed that lobule enrichment for this cluster was strongly consistent with Allen Brain Atlas expression staining (Extended Data Fig. 3c).

Continuity of gene expression

To characterize molecular variation across cell types, we attempted to quantify the continuity of scaled gene expression across a given cell type pair, ordered by pseudotime rank (calculated using Monocle2). For each gene, we fit a logistic curve to the scaled gene expression values and calculated the maximum slope (m) of the resulting curve, after normalizing for both the number of cells and dynamic range of the logistic fit. To limit computational complexity, we downsampled cell type pairs to 5,000 total nuclei.

We fit curves and computed m values for the most significantly differentially expressed genes across five cell type pairs (Fig. 3b). Differentially expressed genes were determined using Seurat’s FindMarkers function. We then plotted the cumulative distribution of m values for the top 200 genes for each cell type pair; genes were selected based on ordering by absolute Spearman correlation between scaled gene expression and pseudotime rank.

Trajectory analysis of peri- and postnatal mouse cerebellum data

After generation of digital gene expression matrices for the peri- and postnatal mouse profiles, we filtered out nuclei with fewer than 500 UMIs. We applied the LIGER workflow (similarly to the adult mouse data analysis), to identify clusters corresponding to major developmental pathways. We then isolated the cluster corresponding to GABAergic progenitors (marked by expression of Tfap2b and other canonical markers). We performed a second iteration of LIGER iNMF and Louvain clustering on this population and generated a UMAP representation. Using this UMAP representation, we calculated pseudotime ordering and a corresponding trajectory graph with Monocle3⁴¹. To identify modules of genes which varied along the computed trajectory, we used the graph_test and find_gene_modules functions from Monocle3.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.