Phenotypic variation of transcriptomic cell types in mouse motor cortex

Nature

No statistical methods were used to predetermine sample size. The experiments were not randomized and investigators were not blinded to allocation during experiments and outcome assessment, unless otherwise stated.

Animals

Experiments on adult male and female mice (n = 266; median age 75 days, interquartile range 64–100, full range 35–245 days, Extended Data Fig. 2a) were performed on wild-type C57Bl/6 (n = 27), Viaat-Cre/Ai9 (vesicular inhibitory amino acid transporter, encoded by the Slc32a1 gene, n = 24), Sst-Cre/Ai9 (somatostatin, n = 75), Vip-Cre/Ai9 (vasoactive intestinal polypeptide, n = 46), Pvalb-Cre/Ai9 (parvalbumin, n = 76), Npy-Cre/Ai9 (neuropeptide Y, n = 2), Vipr2-Cre/Ai9 (vasoactive intestinal peptide receptor 2, n = 7), Scl17a8-Cre/Ai9 (VGLUT3, vesicular glutamate transporter 3, n = 6), Gnb4-Cre/Ai9 (n = 1), and Slc17a8-iCre/Ai9 (n = 2) mice. Numbers above refer to mice from which sequencing data were successfully obtained. Several more animals were used for measuring layer boundaries and follow-up experiments at physiological temperature (see below). Mice were co-housed with littermates (2–5 per cage) in a controlled environment at 22–24 °C and 30–70% humidity. Mice were maintained with unrestricted access to food and water on a 12-h light/dark cycle. Procedures for mouse maintenance and surgeries were performed according to protocols approved by the Institutional Animal Care and Use Committee (IACUC) of Baylor College of Medicine.

The Viaat-Cre line was generously donated by Huda Zoghbi (Baylor College of Medicine), the Slc17a8-iCre line by Rebecca Seal (University of Pittsburg). The Gnb4-Cre line was  from the Allen Institute for Brain Science. The other Cre and reporter lines were purchased from the Jackson Laboratory: Sst-Cre (stock no. 013044), Vip-Cre (stock no. 010908), Pvalb-Cre (stock no. 008069), Vipr2-Cre (stock no. 031332), Slc17a8-Cre (stock no. 028534), Npy-Cre (stock no. 027851), Ai9 reporter (stock no. 007909).

We were unable to find any labelled cells in MOp in the Gnb4-Cre mice: all labelled cells were far outside of MOp and close to the claustrum46. For this reason, the data set does not include any Gnb4-positive cells.

Slice preparation

The MOp brain slices were obtained following previously described protocols5,28. In brief, the animals were deeply anaesthetized using 3% isoflurane and decapitated. The brain was rapidly removed and collected into cold (0–4 °C) oxygenated NMDG (N-methyl-d-glucamine) solution containing 93 mM NMDG, 93 mM HCl, 2.5 mM KCl, 1.2 mM NaH2PO4, 30 mM NaHCO3, 20 mM HEPES, 25 mM glucose, 5 mM sodium ascorbate, 2 mM thiourea, 3 mM sodium pyruvate, 10 mM MgSO4 and 0.5 mM CaCl2, pH 7.35 (all from Sigma-Aldrich). We cut 300-μm-thick coronal slices using a Leica VT1200 microtome following coordinates provided in the Allen Brain Atlas for adult mouse (http://atlas.brain-map.org). The slices were subsequently incubated at 34.0 ± 0.5 °C in oxygenated NMDG solution for 10–15 min before being transferred to the artificial cerebrospinal fluid (ACSF) solution containing: 125 mM NaCl, 2.5 mM KCl, 1.25 mM NaH2PO4, 25 mM NaHCO3, 1 mM MgCl2, 11.1 mM glucose and 2 mM CaCl2, pH 7.4 (all from Sigma-Aldrich) for about 1 h. The slices were allowed to recover in ACSF equilibrated with CO2/O2 gas mixture (5% CO2, 95% O2), at room temperature (approximately 25 °C) for 1 h before experiments. During the recordings, slices were submerged in a customized chamber continuously perfused with oxygenated physiological solution. Recorded cells were generally located 15–60 μm deep under the slice surface.

Patch-seq recording procedures

In order to simultaneously obtain electrophysiological, morphological and transcriptomic data from the same neurons, we applied our recently developed Patch-seq protocol17, with some modifications. In particular, changes were made to the internal solution to optimize its osmolarity in order to improve staining quality. RNase-free intracellular solution was prepared as follows: we dissolved 111 mM potassium gluconate, 4 mM KCl, 10 mM HEPES and 0.2 mM EGTA in RNase-free water in a 125-ml Erlenmeyer flask. We then covered the solution with aluminium foil and autoclaved it. After the solution was cooled down to room temperature, we added 4 mM MgATP, 0.3 mM Na3GTP, 5 mM sodium phosphocreatine, and 13.4 mM biocytin (all from Sigma-Aldrich). The pH was adjusted to 7.25 with RNase-free 0.5 M KOH using a dedicated pH meter (cleaned with RNase Zap and RNase-free water before each use). RNase-free water was then added to the solution in order to obtain the desired volume. After carefully checking its osmolarity (approximately 235–240 mOSM) the solution was stored at −20 °C and used for no longer than 3 weeks.

Before each experiment, we combined 494 μl internal solution with 6 μl recombinant RNase inhibitor (1 U/μl, Takara) to increase RNA yield. The addition of the inhibitor resulted in an increase in osmolarity to the desired value of 315–320 mOSM without a further dilution17. The osmolarity of the ACSF was monitored before each experiment and adjusted to be 18–20 mOSM lower than the internal solution. In particular, when the ACSF osmolarity was too low, we added a small amount of sucrose to ACSF to increase its osmolarity and bring it to the desired range. This osmolarity difference between ACSF and the internal solution is important to obtain slight swelling of the cell during the recording session, which improves the diffusion of biocytin in the neuronal processes. All glassware, spatulas, stir bars, counters, and anything else that may come into contact with the reagents or solution were cleaned thoroughly with RNase Zap before use.

Recording pipettes (B200-116-10; Sutter Instrument) of ~3–7 MΩ resistance were filled with 0.1–0.3 μl RNase-free intracellular solution. The size of the pipette tip was chosen according to the target neuron size: 3–4-MΩ pipettes were used to record large neurons (for example, L5 ET excitatory neurons) and 6–7-MΩ pipettes were used to record small cells such as L1 or Vip interneurons.

The PatchMaster software (HEKA Elektronik) and custom Matlab scripts were used to operate the Quadro EPC 10 amplifiers and to perform online and offline data analysis. We used the following quality control criteria: (1) seal resistance value >1 GΩ before achieving whole-cell configuration; (2) access resistance <30 MΩ. Each neuron was injected with 600-ms-long current pulses starting from −200 pA and up to 1,380 pA with 20-pA increment steps (in some cases stimulation was stopped before reaching 1,380 pA). There were 1.3- or 1.4-s intervals between successive current pulses, depending on the used setup. For most neurons, the stimulation was then repeated multiple times from the beginning. Electrophysiological traces used for the analysis were acquired between 3 and 15 min after achieving the whole-cell configuration. Recordings were performed at room temperature (25 °C), as opposed to physiological temperature (34 °C), in order to keep the cells alive for longer. We performed control experiments at physiological temperature as well (see below).

Typically, excitatory neurons were recorded for 5–20 min while interneurons were recorded for 20–50 min in order to allow biocytin to diffuse into distal axonal segments. During the recording, the access resistance was checked every three minutes in order to maintain a stable seal that would ensure successful biocytin diffusion. The resulting cDNA yield was not correlated with the hold time (Spearman correlation −0.01).

Experiments at physiological temperature

A subset of electrophysiological recordings was performed at 34 °C in the presence of fast glutamatergic and GABAergic synaptic transmission blockers, 1 mM kynurenic acid (Sigma-Aldrich) and 0.1 mM picrotoxin (Tocris), respectively. The temperature was maintained stable, and constantly monitored using the temperature controller TC07 (Luigs and Neumann). In this set of experiments, the morphologies were not recovered and multiple neurons were recorded in each slice. The soma depth and the slice thickness were measured before each recording using Linlab2 software (Scientifica). Intrinsic electrophysiological recordings were obtained using the same stimulation paradigm as described above.

In these experiments, we targeted L5 Sst and excitatory neurons (Extended Data Fig. 6). We sequenced in total 185 neurons, obtained from 8 adult mice (7 Sst-Cre/Ai9 and 1 Pvalb-Cre/Ai9), of which 177 neurons passed the transcriptomic quality control and got a t-type assignment (see below). One hundred and ten cells mapped to the Sst subclass, 43 to IT, 12 to ET, 10 to Pvalb, and 2 to NP. 175 cells were assigned to L5 in the post hoc analysis (see below). We obtained high-quality electrophysiological recordings and extracted electrophysiological features of 184 cells.

RNA sequencing of patched cells

At the end of the recording session, cell contents were aspirated into the glass pipette by applying a gentle negative pressure (0.7–1.5 pounds per square inch) for 1–5 min until the size of the cell body was visibly reduced. In most cases, the cell nucleus was visibly attached to the pipette tip and extracted from the cell body. We avoided complete nucleus aspiration, because it can lead to the collapse of the soma structure and of the nearby neurites, resulting in lower staining quality and stronger background staining. During the aspiration process, the cell body structure and access resistance were constantly monitored. Special care was taken to ensure that the seal between the pipette and the cell membrane remained intact to reduce contamination from the extracellular environment. After aspiration, the contents of the pipette were immediately ejected into a 0.2-ml PCR tube containing 4 μl lysis buffer (with ERCC spike-ins), and RNA was subsequently converted into cDNA using a Smart-seq2-based protocol21 as described previously17. The resulting cDNA libraries were screened using an Agilent Bioanalyzer 2100. Samples containing less than around 1 ng total cDNA (in the 15 μl final volume) or with an average size less than 1,500 bp were typically not sequenced (with some occasional exceptions). The cDNA libraries were then frozen and sent for sequencing in 12 separate batches.

The cDNA libraries derived from each neuron were purified and 0.2 ng of the purified cDNA was tagmented using the Illumina Nextera XT Library Preparation with one-fifth of the volumes stated in the manufacturer’s recommendation. Custom 8-bp index primers were used at a final concentration of 0.1 μM. The resulting cDNA library of each batch was sequenced on an Illumina NextSeq500 instrument with a sequencing setup of 75-bp single-end reads and 8-bp index reads. The investigators were blinded to the cell type of each sample during library construction and sequencing.

The sequencing data were processed using the zUMIs 2.5.6b pipeline with default settings47. Sequencing reads were aligned to the mm10 mouse reference genome using STAR version 2.5.4b48 and transcript assignment performed with Gencode transcript annotations, version M23. A substantial portion of the RNA extracted from the neurons was nascent and contained intronic sequences. To accommodate this, gene expression counts were separately calculated using reads mapping to annotated intronic and exonic regions. We detected 42,466 genes, including pseudogenes and annotated non-coding segments, in at least one cell. The resulting exonic and intronic read count data were used for all transcriptomic analyses presented here. To quantify gene expression, we typically normalized exon and intron counts by exonic and intronic gene lengths in kilobases and added normalized counts together to obtain normalized exonic + intronic expression levels. See below for more details. Throughout the manuscript, ‘detected gene’ refers to a gene with a non-zero exonic or intronic count.

Biocytin staining and morphological reconstructions

Morphological recovery was carried out as previously described5,17,28. In brief, after the recordings, the slices were immersed in freshly prepared 2.5% glutaraldehyde, 4% paraformaldehyde solution in 0.1 M PBS at 4 °C for at least 48 h. The slices were subsequently processed with the avidin-biotin-peroxidase method to reveal the morphology of the neurons. As previously described, we took several steps to improve the staining quality of the fine axonal branches of interneurons5,17. First, we used a high biocytin concentration (0.5 g/100 ml). Second, we incubated with avidin–biotin complex and detergents at a high concentration (Triton X-100, 5%) for at least 24 h before staining with 3,3′-diaminobenzidine (DAB).

Recovered cells were manually reconstructed using a 100 × oil-immersion lens and a camera lucida system (MicroBrightField). We aimed to reconstruct all cells that had staining of sufficient quality (axons and dendrites for the inhibitory neurons; only dendrites for the excitatory neurons), and obtained 646 reconstructions in total. In addition, we reconstructed the dendrites of 30 neurons from the Vip and Scng subclasses that lacked sufficient axonal staining. Vip neurons are traditionally classified on the basis of dendritic morphology, so these reconstructions can inform t-type characterizations. These additional 30 reconstructions are shown, together with the main 646 reconstructions, in Supplementary File 1.

Forty-five sequenced cells were mistakenly recorded using a solution with a much smaller concentration of biocytin, and their morphologies could not be recovered. We made sure that the measured electrophysiological properties of these cells were not systematically different from those of the the other sequenced cells.

Inevitably, neuronal structures can be severed as a result of the slicing procedure. We took special care to exclude reconstructions of all neurons that showed any signs of damage, lack of contrast, or poor overall staining. Consistently with previous studies, tissue shrinkage due to the fixation and staining procedures was about 10–20%5,28,49. This shrinkage was not compensated for in our analysis.

Cortical thickness normalization and layer assignment

Nissl-stained slices (n = 15 from two wild-type adult mice) were used to measure normalized layer boundaries in MOp. The Nissl staining protocol was adapted from ref. 50. In brief, brain slices were mounted on slides and allowed to dry. The sections were then demyelinated, stained with 0.1% cresyl violet-acetate (C5042, Sigma) for 30 min at 60 °C and further destained. The sections were then coverslipped in Cytoseal 60 (Richard Allan Scientific). For each slice we measured total thickness from pia to white matter and the depths of the three between-layer boundaries (L1 to L2/3, L2/3 to L5, L5 to L6), based on the cortical cytoarchitecture, using a Neurolucida system with 10 × or 20 × magnification. All measurements were normalized by the respective slice thickness, and the averages over all n = 15 slices were used as the normalized layer boundaries (Extended Data Fig. 2b).

For the Patch-seq neurons, we measured soma depth and the cortical thickness of the slice using a Neurolucida system. We took their ratio as the normalized soma depth, and assigned each neuron to a layer (L1, L2/3, L5, or L6) based on the Nissl-determined layer boundaries (Extended Data Fig. 2b). We obtained soma depth information for 1,284 neurons out of 1,329 (45 neurons were mistakenly recorded using a solution with insufficient biocytin concentration, and we could measure soma depths for only 2 of those; for 2 other neurons the measurements could not be carried out because the slices were lost). For the 45 neurons with missing soma depth measurements, we used the layer targeted during the recording for all layer-based analyses and visualizations (marker shapes in Figs. 1c–e, 3a–c, layer-restricted analysis in Fig. 4, Extended Data Fig. 8).

All reconstructed morphologies were normalized by the cortical thickness of the respective slice to make it possible to display several morphologies next to each other, as in Extended Data Fig. 3.

t-Type assignment

The t-type assignment procedure was done in two rounds. The first round was for quality control and initial assignment to one of the three large transcriptomic groups (CGE-derived interneurons, MGE-derived interneurons, and excitatory neurons) that are perfectly separated from each other with no transcriptomically intermediate cells4. The second round was done to assign the cells to specific t-types.

In the first round, we mapped each Patch-seq cell to a large annotated Smart-seq2 reference data set from adult mouse cortex4, using a procedure similar to the one described in ref. 28. Specifically, using the exon count matrix of the reference data set, we selected the 3,000 most variable genes (see below). We then normalized all exon counts by exonic gene lengths in kilobases, all intron counts by intronic gene lengths in kilobases (plus 10−6, to avoid division by zero) and added normalized counts together to obtain normalized exonic + intronic expression levels. We log-transformed these values using log2(x + 1) transformation and averaged the log-transformed values across all cells in each of the 133 t-types, to obtain reference transcriptomic profiles of each t-type (133 × 3,000 matrix). Out of these 3,000 genes, 2,666 were present in the genome annotation that we used and were detected in our data set. We applied the same normalization and log-transformation procedure to the exonic and intronic read counts of our cells, and for each cell computed Pearson correlation across the 2,666 genes with each of the 133 t-types. Each cell was assigned to the t-type to which it had the highest correlation (Extended Data Fig. 1d).

Cells meeting any of the following exclusion criteria were declared low quality and did not get a t-type assignment (Extended Data Fig. 2e): cells with the highest correlation below 0.4 (78 cells); cells that would be assigned to non-neural t-types, presumably owing to RNA contamination51 (14 cells; see also Extended Data Fig. 2j–n); cells with the highest correlation less than 0.02 above the maximal correlation in one of the other two large transcriptomic groups (5 cells). The remaining 1,232 cells passed quality control and entered the second round.

In the second round, cells were independently mapped to the seven transcriptomic data sets obtained from mouse MOp20. The mapping was done only to the t-types from the transcriptomic group identified in the first round, using the 500 most variable genes in that data set for that transcriptomic group (so using 7 × 3 = 21 sets of 500 most variable genes). Gene selection was performed as described below, and t-type assignment was done exactly as described above. Across the 21 reference subsets, 421–494 most variable genes were present in our data set, and were used for the t-type assignment (Extended Data Fig. 1e). When mapping to the Smart-seq2 reference data sets, we used normalized intronic and exonic reference counts, as above. When mapping to the UMI-based reference data sets, we used the unique molecular identifier (UMI) counts directly, without gene length normalization.

We used bootstrapping over genes to assess the confidence of each t-type assignment. For each cell and for each of the seven reference data sets, we repeatedly selected a bootstrap sample of genes (that is, the same number of genes, selected randomly with repetitions) and repeated the mapping. This was done 100 times and the fraction of times the cell mapped to each t-type was taken as the t-type assignment confidence for that t-type (Extended Data Fig. 1f). The confidences obtained with seven reference data sets agreed well with each other (Extended Data Fig. 2i) and were averaged to obtain the consensus confidence. Finally, the cell was assigned to the t-type with the highest consensus confidence.

Four cells were assigned to an excitatory t-type, despite having clearly inhibitory firing, morphology, and/or soma depth location (such as L1). The most likely cause of this was RNA contamination from excitatory cells, which are much more abundant than inhibitory cells in the mouse cortex (Extended Data Fig. 2). These four cells were excluded from all analyses and visualizations (as if they did not pass the transcriptomic quality control). In addition, one cell was probably located outside MOp, based on the slice anatomy, and was excluded as well. The final number of cells with t-type assignment was 1,227.

Selection of most variable genes

Several steps of our analysis required selecting a set of the most variable genes in a given transcriptomic data set. We always selected a fixed predefined number of genes (such as 500, 1,000, or 3,000).

To select the most variable genes, we found genes that had, at the same time, high non-zero expression and a high probability of near-zero expression52. Our procedure is described in more detail elsewhere23. Specifically, we excluded all genes that had counts of at least cmin (for Patch-seq and Smart-seq2: cmin = 32; for 10x: cmin = 0) in fewer than 10 cells. For each remaining gene we computed the mean log2 count across all counts that were larger than cmin (non-zero expression, μ) and the fraction of counts that were smaller than or equal to cmin (probability of near-zero expression, τ). Across genes, there was a clear inverse relationship between μ and τ, that roughly followed the exponential law: τ ≈ exp(−1.5 × μa) for some horizontal offset a. Using a binary search, we found a value b of this offset that yielded the desired number of genes with τ > exp(−1.5 × μb) + 0.002.

For Smart-seq2 and Patch-seq data sets, we used only exonic counts to perform gene selection.

t-SNE visualization of the transcriptomic data

t-SNE embeddings22 of the three subsets of the single-cell 10x v2 data set20 (Fig. 1c–e) were constructed using the same 500 most variable genes that were used for t-type assignment (see above). The UMI counts were normalized by each cell’s sequencing depth (sum of counts), multiplied by the median sequencing depth across all cells, log2(x + 1)-transformed, and reduced to 50 principal components. The resulting n × 50 matrix was used as input to t-SNE. We used FIt-SNE 1.2.153 with default parameters (including learning rate n/12 and scaled principal component analysis (PCA) initialization23). Perplexity was left at the default value of 30 for both inhibitory subsets and increased to 100 for the excitatory subset.

To position Patch-seq cells on a reference t-SNE embedding, we used a published procedure23. In brief, each cell was positioned at the median embedding location of its ten nearest neighbours, based on Pearson correlation distance in the high-dimensional space. As above, we used the sum of the normalized exonic and intronic counts for Patch-seq cells, and raw UMI counts for the reference cells. All values were log2(x + 1)-transformed and correlations were computed across the same genes that were used for t-type assignments (see above).

Extraction of electrophysiological features

Twenty-nine electrophysiological properties of the neurons were automatically extracted based on the raw membrane voltage traces (Extended Data Fig. 4) using Python scripts from the Allen Software Development Kit (SDK) (https://github.com/AllenInstitute/AllenSDK) with some modifications to account for our experimental paradigm (https://github.com/berenslab/EphysExtraction).

For each hyperpolarizing current injection, the resting membrane potential was computed as the mean membrane voltage during 100 ms before stimulation onset and the input resistance as the difference between the steady state voltage and the resting membrane potential, divided by the injected current value (we took the average voltage of the last 100 ms before stimulus offset as steady state). The median of these values over all hyperpolarizing traces was taken as the final resting membrane potential and input resistance, respectively.

To estimate the rheobase (the minimum current needed to elicit any spikes), we used robust regression (random sample consensus algorithm, as implemented in sklearn.linear_model.RANSACRegressor) of the spiking frequency onto the injected current using the five lowest depolarizing currents with non-zero spike count (if there were fewer than five, we used those available). The point at which the regression line crossed the x-axis gave the rheobase estimate (Extended Data Fig. 4). We restricted it to be between the highest injected current that elicited no spikes and the lowest injected current that elicited at least one spike. If the regression line crossed the x-axis outside this interval, the first current step that elicited at least one spike was used.

The action potential (AP) threshold, AP amplitude, AP width, afterhyperpolarization (AHP), afterdepolarization (ADP), the first AP latency, and the upstroke-to-downstroke ratio (UDR) were computed as illustrated in Extended Data Fig. 4, using the first AP fired by the neuron. AP width was computed at the AP half-height. UDR refers to the ratio of the maximal membrane voltage derivative during the AP upstroke to the maximal absolute value of the membrane voltage derivative during the AP downstroke. We also computed the first AP latency at 20 pA current above the smallest current stimulation value that elicited a spike.

The interspike interval (ISI) adaptation index for each trace was defined as the ratio of the second ISI to the first one. The ISI average adaptation index was defined as the mean of ISI ratios corresponding to all consecutive pairs of ISIs in that trace. For both quantities we took the median over the five lowest depolarizing currents that elicited at least three spikes (if fewer than five were available, we used all of them). AP amplitude adaptation index and AP amplitude average adaptation index were defined analogously to the two ISI adaptation indices, but using the ratios of consecutive AP amplitudes (and using the median over the five lowest depolarizing currents that elicited at least two spikes).

The maximum number of APs refers to the number of APs emitted during the 600-ms stimulation window of the highest firing trace. The spike frequency adaptation (SFA) denotes the ratio of the number of APs in the second half of the stimulation window to the number of APs in the first half of the stimulation window of the highest firing trace. If the highest firing trace had fewer than five APs, SFA was not defined. Here and below the highest firing trace corresponds to the first depolarizing current step that showed the maximum number of APs during the current stimulation window (after excluding all stimulation currents for which at least one AP was observed in 100 ms before or in 200 ms after the stimulation window; see below).

The membrane time constant (τ) was computed as the time constant of the exponential fit to the membrane voltage from the stimulation onset to the first local minimum (we took the median over all hyperpolarizing traces). Three further features described the sag of the first (the lowest) hyperpolarization trace. The sag ratio was defined as the difference between the sag trough voltage (average voltage in a 5-ms window around the sag trough) and the resting membrane potential, divided by the steady state membrane voltage difference from the resting membrane potential. The sag time was defined as the time period between the first and the second moments at which the membrane voltage crossed the steady-state value after the stimulation onset. The sag area refers to the absolute value of the integral of the membrane voltage minus the steady-state voltage during the sag time period (Extended Data Fig. 4). If the sag trough voltage and the steady-state voltage differed by less than 4 mV, the sag time and sag area were set to zero.

The rebound was defined as the voltage difference between the resting membrane potential and the average voltage over 150 ms (or whatever time remained until 300 ms after the stimulation offset) after rebound onset, which we identified as the time point after stimulation offset at which the membrane voltage reached the value of the resting membrane potential. If the membrane voltage never reached the resting membrane potential during the 300 ms after the stimulation offset, the rebound was set to zero. The rebound number of APs was defined as the number of APs emitted during the same period of time. Both rebound features were computed using the lowest hyperpolarization trace.

The ISI coefficient of variation (CV) refers to the standard deviation divided by the mean of all ISIs in the highest firing trace. Note that a Poisson firing neuron would have ISI CV equal to one. The ISI Fano factor refers to the variance divided by the mean of all ISIs in the highest firing rate. The AP CV and AP Fano factor refer to the CV and the Fano factor of the AP amplitudes in the highest firing trace, respectively.

The burstiness was defined as the difference between the inverse of the smallest ISI within a detected burst and the inverse of the smallest ISI outside bursts, divided by their sum. We took the median over the first five depolarizing traces. We relied on the Allen SDK code to detect the bursts. In brief, within that code a burst onset was identified whenever a ‘detour’ ISI was followed by a ‘direct’ ISI. Detour ISIs are ISIs with a non-zero ADP or a drop of at least 0.5 mV of the membrane voltage after the first AP terminates and before the next one is elicited. Direct ISIs are ISIs with no ADP and no such drop of membrane voltage before the second AP. A burst offset was identified whenever a direct ISI was followed by a detour ISI. Additionally, bursts were required to contain no ‘pauselike’ ISIs, defined as unusually long ISIs for that trace (see Allen SDK for the implementation details).

Some neurons (in particular neurogliaform cells) started to emit APs before and after the current stimulation window, after the stimulation currents exceeded a certain amount. To quantify this effect, we defined wildness as the difference in the number of APs between the highest firing trace (possibly showing APs before or after the stimulation window) and the highest firing trace as defined above (without any APs outside the stimulation window). For most neurons, wildness was equal to zero.

For all statistical analysis we used 17 features out of the extracted 29, excluding features that were equal to zero for many cells (afterdepolarization, burstiness, rebound number of APs, sag area, sag time, wildness), two Fano factor features that were highly correlated with the corresponding coefficient of variation features (AP Fano factor, ISI Fano factor) and another measure of latency that was highly correlated with the latency itself, features that had very skewed distributions (AP amplitude average adaptation index, ISI average adaptation index), and features that were undefined for some of the cells (spike frequency adaptation). Four features were log-transformed to make their distribution more Gaussian-like: AP coefficient of variation, ISI coefficient of variation, ISI adaptation index, and latency.

Extraction of morphological features

Reconstructed morphologies were converted into the SWC format using NLMorphologyConverter 0.9.0 (http://neuronland.org) and further analysed using MorphoPy (https://github.com/berenslab/MorphoPy, version 0.6)54. Each cell was soma-centred in the x (slice width) and y (slice depth) dimensions, and aligned to pia in the z (cortical depth) dimension so that z = 0 corresponded to pia. All neurites were smoothed in the slice depth dimension (y) using a Savitzky–Golay filter of order 3 and window length 21, after resampling points to have maximally 1 μm spacing. For further analysis we computed two different feature representations of each cell: the normalized z-profile and a set of morphometric statistics24,28,55.

To compute the normalized z-profile, we divided all the coordinates of the neuronal point cloud by the thickness of the respective cortical slice, so that z = 1 corresponded to the white matter border. We projected this point cloud onto the z-axis and binned it into 20 equal-sized bins spanning [0, 1]. The resulting histogram describes a neuron’s normalized depth profile perpendicular to the pia. For the purposes of downstream analysis, we treated this as a set of 20 features. The z-profiles were separately computed for axons and dendrites.

Morphometric statistics were separately computed for the dendritic and axonal neurites to quantify their arborization shape and branching patterns. For the excitatory neurons, several additional morphometric statistics were computed for the apical dendrites, where apical dendrite was operationally defined as the dendrite with the longest total path length. We further used two ‘somatic’ features: normalized soma depth and soma radius. We did not use any features measuring morphological properties in the slice depth (y) direction because of possible slice shrinkage artefacts. We did not use any axonal features for the excitatory cells because only a small part of the axon could typically be reconstructed. For the inhibitory cells, where dendrite and axon could both be fully recovered, we included some measures of dendritic and axonal overlap. The full list of morphometric statistics is given in Supplementary File 3.

We extracted a set of 75 features, of which 40 were defined for excitatory neurons and 62 for inhibitory neurons, and processed the data for excitatory and inhibitory neurons separately. In each case, we excluded features with coefficient of variation below 0.25 (among the features with only positive values). This procedure excluded five features for the excitatory and nine features for the inhibitory cells. The distributions of the remaining features were visually checked for outliers and for meaningful variation between transcriptomic types, leading to a further exclusion of three features for the inhibitory cells. The full list of excluded features is given in Supplementary File 3. The resulting set of morphometric statistics used for further analysis consisted of 35 features defined for the excitatory neurons and 50 features defined for the inhibitory neurons.

Reduced-rank regression

For the RRR analysis32 we used 17 electrophyiological features and all 1,219 cells for which values for all 17 features and a t-type assignment could be computed. Electrophysiological features were standardized. Exon counts and intron counts were normalized by the exon/intron gene lengths as described above, summed together, converted to CPM, log2(x + 1)-transformed, and then standardized. We selected the 1,000 most variable genes (using raw exonic counts) and used only those for the RRR analysis.

In brief, RRR finds a linear mapping of gene expression levels to a low-dimensional latent representation, from which the electrophysiological features are then predicted with another linear transformation (for mathematical details, see ref. 32). The model uses sparsity constraints in the form of elastic net penalty to select only a small number of genes. For Fig. 2 we used a model with rank r = 5, zero ridge penalty (α = 1), and lasso penalty tuned to yield a selection of 25 genes (λ = 0.5). Cross-validation (Extended Data Fig. 5) was done using 10 folds, elastic net α-values 0.5, 0.75, and 1.0, and λ-values from 0.2 to 6.0.

The plots shown in Fig. 2a, b are called bibiplots because they combine two biplots: the left biplot shows a mapping of gene expression levels onto the two latent dimensions; the right biplot shows the same mapping of electrophysiological features. To illustrate the meaning of the latent dimensions, each biplot combines the resulting scatter plots with lines showing how original features are related to the latent dimensions. Specifically, we computed the correlations of individual genes or electrophysiological properties with the latent dimensions and visualized these correlations as lines on the biplot. The circle shows the maximal possible correlation; only lines longer than 0.4 times the circle radius are shown in Fig. 2. Label positions were automatically adjusted by simulating repulsive forces between all overlapping pairs of labels, until there was no overlap.

For the model based on ion channel genes, we obtained the list of 328 ion channel genes from https://www.genenames.org/data/genegroup/#!/group/177and used all 307 of them that had non-zero expression in at least 10 of our cells. We used rank r = 5, α = 1, and λ tuned to yield 25 genes (λ = 0.303), as above.

t-SNE visualization of the morpho-electric phenotypes

For the t-SNE visualization22 of the electrophysiological phenotypes, we used 17 features as described above and all n = 1,320 cells that had values for all 17 features. All features were standardized across this set of cells and transformed with PCA into a set of 17 PCs. We scaled the PCs by the standard deviation of PC1. We used the t-SNE implementation from scikit-learn Python library with the default perplexity (30), early exaggeration 4 (the default value 12 can be too large for small data sets), and scaled PCA initialization23. Fig. 3a shows n = 1,219 cells that had a t-type assignment.

For the t-SNE visualization of the morphological phenotypes, we combined morphometric statistics with the normalized z-profiles. The pre-processing, including PCA, was done separately for the excitatory and inhibitory neurons because they used different sets of morphometric statistics (see above). Only neurons with assigned t-types were used for this analysis. Two inhibitory neurons were left out because some of the morphometric statistics could not be extracted owing to insufficient dendritic recovery; this left 367 inhibitory neurons (with 50 morphometric features) and 269 excitatory neurons (with 35 morphometric features). All features were standardized and each set was reduced to 20 PCs. We scaled the PCs by the standard deviation of the respective PC1, to make the inhibitory and the excitatory PCs have comparable variances.

We used dendritic z-profiles for the excitatory neurons and axonal z-profiles for the inhibitory neurons. We reduced each set to five PCs, discarded PC1 (it was strongly correlated with the normalized soma depth and made the resulting embedding strongly influenced by the soma depth), and scaled the PCs by the standard deviation of the respective PC2. We stacked the 20 scaled morphometric PCs and the 4 scaled z-profile PCs to get a combined 24-dimensional representation, separately for the excitatory and for the inhibitory neurons. We then combined these representations into one block-diagonal 48-dimensional matrix. This procedure makes the excitatory and the inhibitory populations both have zero mean. To prevent overlap between these two populations, we added a small constant value of 0.25 to the excitatory block-diagonal block, leading to the strong excitatory–inhibitory separation in Fig. 3b. The t-SNE was performed exactly as described above.

For the t-SNE visualization of the morpho-electrical landscape, we stacked together the 48-dimensional morphological representation and the 16-dimensional electrophysiological representation obtained above, using only cells that had all morphological and all electrophysiologcal features (n = 628). We multiplied the electrophysiological block by √2 to put its total variance on a similar scale (it only consisted of one set of scaled PCs, whereas the morphological representation consisted of two sets of scaled PCs: morphometrics and z-profiles). The resulting 64-dimensional morpho-electrical representation was used for t-SNE, exactly as described above.

kNN classification of transcriptomic families

To classify neurons into transcriptomic families on the basis of electrophysiological, morphological, or combined features (Figs. 3d, 5a, Extended Data Fig. 8a), we used a kNN classifier with k = 10 and Euclidean distance metric (taking the majority family among the k nearest neighbours). This is effectively a leave-one-out cross-validation procedure. For each data modality we took the exact same data representation that was used for computing t-SNE embeddings (Fig. 3a–c; see above). Note that the t-SNE algorithm is also based on nearest neighbours and makes all close neighbours attract each other in the embedding. We chose the kNN classifier as a simple but versatile non-parametric classifier that is directly related to the t-SNE embeddings. We did not use the Sncg and NP families owing to insufficient coverage in our data set (Fig. 1).

Fig. 3d shows the fraction of cells from each family that was classified into each family. Fig. 5a and Extended Data Fig. 8a show fractions of cells from each t-type that were classified into each family. For morphological and combined features, Extended Data Fig. 8a shows fractions of cells from the majority layer of each t-type. For example, the Pvalb Reln type occurred most often in L5, so only cells from that layer were taken for that type. Only t-types with at least ten cells (or at least ten layer-restricted cells) are shown.

Within-family analysis

To study the relationship between transcriptomic and electrophyiological distances between pairs of t-types (Fig. 4c, d, Extended Data Figs. 6, 7), we took all t-types with five or more cells assigned to them (for Extended Data Fig. 7a: with ten or more). For each pair of t-types, transcriptomic distance was computed as the Pearson correlation between the average log2(x + 1)-transformed UMI counts in the single-cell 10x v2 data20. The 1,000 most variable genes across all neural types were used for Fig. 4c and Extended Data Fig. 7a, b, and the 500 most variable genes across the respective transcriptomic group (see above) were used for Fig. 4d and Extended Data Figs. 6i, j and 7c–n. Electrophysiological distance was computed as the Euclidean distance between the average feature vectors. Fig. 4d used the soma depth distance, computed as the absolute value of the difference between the average normalized soma depths.

T-type variability analysis

The normalized total variance in Fig. 5b and Extended Data Fig. 8b was computed as follows. For each modality, we took the exact same data representation that was used for computing t-SNE embeddings (Fig. 3a–c; see above). For each t-type (or layer-restricted t-type; see above), we took the sum of its variances in all dimensions as the total variance and divided by the sum of variances in all dimensions across the whole data set:

$$frac{{sum }_{j}frac{1}{|T|}{sum }_{iin T}{left({X}_{ij}-frac{1}{|T|}{sum }_{iin T}{X}_{ij}right)}^{2}}{{sum }_{j}frac{1}{n}{sum }_{i}{left({X}_{ij}-frac{1}{n}{sum }_{i}{X}_{ij}right)}^{2}},$$

where Xij is a value of feature j of cell i, n is the total number of cells, and T is the set of cell numbers belonging to the given t-type. The value 0 indicates that all cells from this t-type have exactly identical features. The value 1 indicates that there is as much variance in this one t-type as in the whole data set. Only t-types with at least ten cells (or at least ten layer-restricted cells) are shown in Fig. 5b and Extended Data Fig. 8b.

To provide a sensible baseline for the range of possible normalized total variances in a population of morpho-electrically homogeneous types, we used a clustering analysis. For the cells of all the K t-types (or layer-restricted t-types) with at least ten cells in a given panel, we used the k-means algorithm to cluster them into K clusters, reasoning that these clusters should be as homogeneous as possible given the variability in our data set. We used the k-means implementation from scikit-learn with default parameters. We then computed the normalized total variance of each cluster as described above. Grey shading in Fig. 5b and Extended Data Fig. 8b shows the interval between the minimum and the maximum cluster variances. Note that the k-means algorithm directly minimizes within-cluster total variances.

We used the entropies of a Leiden clustering35 as an alternative way to approach the same question. For each modality, using the exact same data representation as above, we constructed its kNN graph with k = 10 and clustered it using the Leiden algorithm as implemented in the Python package leidenalg with RBConfigurationVertexPartition quality function and resolution parameter manually tuned to produce roughly the same number of clusters for each modality as in ref. 24. (Extended Data Fig. 8). For each t-type (or layer-restricted t-type), we then measured the entropy of the distribution of electrophysiological or morphological cluster IDs, after randomly subsampling the t-type to ten cells. Subsampling was done to eliminate a possible bias due to the t-type abundance. The whole procedure was repeated 100 times with different random seeds for the Leiden clustering and for the subsampling.

Reporting summary

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

Products You May Like

Articles You May Like

Goldman questions the future and purpose of OPEC+ ahead of a critical meeting
A Sweet Thanksgiving for Our Galaxy
Restaurant tech start-up Toast soars to $8 billion valuation seven months after cutting half its staff
Rare Plant Rises from the Dead in Norfolk Pond
In call with troops, Trump basks in Space Force achievement

Leave a Reply

Your email address will not be published. Required fields are marked *