Integration of electrophysiological recordings with single-cell RNA-seq data identifies novel neuronal subtypes

Abstract

Traditionally, neuroscientists have defined the identity of neurons by the cells’ location, morphology, connectivity and excitability. However, the direct relationship between these parameters and molecular phenotypes establishing cellular heterogeneity has remained largely unexplored. Here, we present a method to obtain full transcriptome data from single neocortical pyramidal cells and interneurons after whole-cell patch-clamp recordings in mouse brain slices. In our approach, termed Patch-seq, a patch-clamp stimulus protocol is followed by aspirating the entire somatic compartment into the recording pipette, amplification by single-cell tagged reverse transcription in a tube reaction and Illumina sequencing. We show that Patch-seq reveals a close link between electrophysiological characteristics, responses to acute chemical challenges and RNA expression of neurotransmitter receptors and channels. Moreover, it distinguishes neuronal subpopulations that correspond to both well-established and, to our knowledge, hitherto undescribed neuronal subtypes. Our findings demonstrate the ability of Patch-seq to precisely map neuronal subtypes and predict their network contributions in the brain The morphology, excitability, connectivity and neurotransmitter utilization of individual neurons underlie the distinct computations each neuronal circuit can perform in the nervous system 1-3 . Thus, the identification of distinct subclasses of neurons remains a key challenge in neuroscience. Neuronal taxonomy based on a combination of developmental, morphological and neurophysiological traits is well accepted, particularly for interneurons in the cerebral cortex 4-6 . These classification systems primarily rely on candidate marker analysis by a mixture of patch-clamp electrophysiology and single-cell semi-quantitative PCR (qPCR) 4,7-10 . More recently, advances in single-cell RNA sequencing (RNA-seq) in the central nervous system 11-13 led to the identification of novel cell-types. Particularly, RNA-seq allowed the molecular classification of neurons in the somatosensory cortex and CA1 subfield of the hippocampus into 47 subtypes, including 16 subclasses of interneurons 11 . Despite pioneering work using microarrays 14-16 , multiplexed qPCR 17,18 and even proof-of-concept RNA-seq on single neurons 19 , no robust method exists to simultaneously investigate the electrophysiology, morphology and transcriptome profiles of the same neuron. Combining patch-clamp electrophysiology and post-hoc morphological reconstructions with the resolution of quantitative RNA-seq in single neurons would present a potentially critical advance for neuronal classification as it can resolve transcriptome-wide variations in gene expression to reveal cell-type-specific determinants of neuronal cytoarchitecture and biophysical properties. Nevertheless, only Qiu et al. 19 attempted RNA-seq on material taken from three neurons through a patch pipette without prior patch-clamp recordings, which yielded RNA-seq data of variable quality. Here, we describe Patch-seq, a method relying on sequencing RNA aspirated from the soma of single patch-clamp-recorded neurons. We validate Patch-seq transcriptomes by aligning them with larger single-cell datasets to achieve high-quality classification, particularly to resolve cortical (inter-)neuron types previously considered homogenous into distinct subtypes. We also show that Patch-seq is compatible with the post-hoc morphological analysis of neurons in optically cleared tissues, and produces a quantitative dataset that simultaneously resolves mRNAs for all known ion channels, receptors and synaptic proteins. Acute pharmacological probing of cortical interneurons established causality between RNA-seq-based predictions and experimentally observed neuronal responses ex vivo. Thus, Patch-seq is suited to discover molecular determinants of neuronal morphology and excitability. Results Data collection and characterization We focused on cholecystokinin (CCK)-containing(+) GABAergic interneurons because their morphological and molecular features are thought to form a quasi-continuum from axon- to dendrite-targeting interneurons in cortical areas 3,5,11 . CCK+ interneurons are the plastic and dynamic gate-keepers of neuronal circuits 20 . Their inactivation likely contributes to anxiety, mood disorders and schizophrenia 21,22 . As reliable histochemical detection of CCK+ interneurons is particularly challenging, we established a dual-labeled CCKBAC/dsRedT3::GAD67gfp/+ mouse reporter 23,24 (Fig. 1a, Supplementary Fig. 1a), and sampled dual DsRed+/GFP+ interneurons in layers (L)1/2 of the somatosensory cortex. We took advantage of moderate CCK expression in cortical pyramidal cells 11 (DsRed+ only) to build a reference database of electrophysiological and molecular features that resolve cortical layer specificity for comparative analysis (Fig. 1a). We first selected ~120 DsRed+ cortical neurons in total for patch-clamp electrophysiology and morphological examination (as confirmed by epifluorescence microscopy prior to recording; Fig. 1a), of which 83 cells proved suitable for combined electrophysiology and RNA-seq analysis. Forty-five were interneurons (inhibitory types, “I-types” ; Fig. 1b-f) and DsRed+/GFP+ with an additional 38 being pyramidal cells (excitatory types, “Exc-types”; Supplementary Fig. 1b-d). We only discarded neurons if their access resistance showed >20% deviation from baseline (Online Methods) or if their electrophysiology parameters did not resemble any I-type cluster (3 cells; data not shown). For each neuron, we recorded their location (cortical layer), and analyzed their passive and active membrane properties through a series of custom-written routines for whole-cell patch-clamp electrophysiology (on a HEKA EPC-10 controlled by PatchMaster 2.80; Supplementary Table 1). We also combined patch-clamp electrophysiology, biocytin filling, chromogenic (DAB) labeling and tissue clearing for the reconstruction of (axonal) morphology in optically cleared tissues 25,26 (Fig 1b-f). Our visualization method is also compatible with light-sheet 27 or two-photon microscopy in whole brain slices (300-350-μm thick), thus significantly reducing processing time and increasing the three-dimensional integrity of neuronal morphology. Functional and morphological analysis of CCK interneurons The most common nomenclature for (inter-)neurons 6 combines key biophysical and morphological features. Based on their pattern of postsynaptic target innervation, interneurons fall into one of three morphologically distinct subclasses: axo-axonic, perisomatic or dendrite-targeting cells 3,5-7,28 . Here, we first used patch-clamp electrophysiology to classify DsRed+/GFP+ interneurons in cortical L1 into 5 subclasses (“I-types 1-5, Fig. 1b-f, Supplementary Table 1 and Online Methods) purely based on their electrophysiological properties: I-type 1 interneurons exhibited action potential (AP) accommodation with a low adaptation ratio. They began to spike at the onset of a rheobasic stimulus without producing a burst, demonstrated the smallest AP amplitude, the largest afterhyperpolarization (AHP) and most hyperpolarized resting membrane potential (Vrest, Fig. 1b). I-type 2 cells had accommodating AP trains without producing burst firing. These cells displayed the highest AP amplitude, the steepest AP upstroke slope and the highest firing frequency of all recorded neurons. They did not demonstrate ’sag’ depolarization indicating the activation of hyperpolarization-activated non-selective cationic currents (Fig. 1c). I-type 3 interneurons produced AP bursts, shallowest AP upstroke and significant accommodation (Fig. 1d). Each AP burst consisted of 3-4 spikes after the first AP on 2-times threshold current. These neurons had a large afterdepolarization (ADP), as well as the largest ’sag’ depolarization. I-type 4 interneurons were accommodating cells that displayed a small ADP before a slow AHP; yet this was insufficient to produce a burst (Fig. 1e). These neurons had the highest input resistance amongst all interneuron subclasses examined. I-type 5 interneurons exhibited irregular spiking, and fired APs without producing a burst (Fig. 1f). These cells had a small amplitude ADP on the falling trajectory of the AHP rise, and had the highest rheobase amongst the interneuron types assessed. For comparison and validation, we examined the 38 pyramidal cells recorded in L2/3, L4 and L5 (“Exc L2/3, L4 and L5”, Supplementary Fig. 1b-d, Supplementary Table 1 and Online Methods). L2/3 and L4 pyramidal cells (DsRed only) were categorized as non-bursting, accommodating and regular spiking (Supplementary Fig. 1b,c). In contrast, L5 pyramidal cells exhibited higher input resistance and pronounced ADP and a 2nd consecutive AP when stimulated at 2x the threshold current (Supplementary Fig. 1d). Next, we used post-hoc morphological reconstructions to identify the axonal nests and dendritic arbors of the 5 putative interneuron subtypes (Fig. 1b-f). Our analysis revealed substantial differences in cell morphology, such as: I-type 1 interneurons had equally-sized processes throughout, thus precluding their unequivocal assignment as axons or dendrites at the light-microscopy level. These processes concentrated in L1 and L2. I-type 1 cells maintained gap junctions with neighboring non-pyramidal cells and even to L2/3 pyramidal cells, as shown by dye-loading experiments using Lucifer yellow (see also Fig. 6a,b). I-type 2 interneurons had somatic diameters at least twice of I-type 1 cells, and resided at the L1-L2 transition. Whereas their dendrites intruded into both L1 and L2, their axons were almost exclusively located in L2. I-type 3 interneurons, similar to I-type 1 cells, also had processes with homogenous diameters yet all were less than half of the process thickness of I-type 1 cells. Their axonal arbor was less branched with the majority of processes coursing in L2. I-type 4 cells had the majority of their elaborate axonal nests concentrated in L2 with some collaterals reaching as far as L4. Their dendrites targeted apically. I-type 5 interneurons branched horizontally in L1 and upper L2. Rarely, we saw dye coupling between I-type 5 interneurons and neighboring pyramidal cells, with less dense network labeling than for I-type 1 cells . Cumulatively, the above correlated electrophysiology and morphological differences identify I-type 1-5 CCK+ interneurons as bona fide interneuron subtypes. Their distinct features suggest a critical diversity of molecular determinants for cell identity, which can be interrogated by single-cell RNA-seq. RNA sequencing of somatic aspirates We next developed a method, called Patch-seq, for performing single-cell RNA sequencing on the same neurons that we characterized by patch-clamp electrophysiology. (Fig. 2, 3a-c and Online Methods). After testing the neurons in a series of current-clamp and voltage-clamp protocols within 20-25 min (see above), their entire somatic compartment was aspirated into the recording pipette. We found, using an iterative approach on ~140 cells (Supplementary Fig. 2), that applying a continuum of positive voltage pulses (to membrane potential +20 mV from holding potential −5 mV, each 5 ms in length with 5 ms intervals) reduced the loss of RNA by most efficiently holding negatively-charged RNA molecules in the pipette solution. Subsequently, the samples (0.8-0.9 μl) were ejected into lysis buffer (0.6 μl; Online Methods). RNA-seq was performed using a modified version of single-cell tagged reverse transcription (STRT-C1) 29,30 designed to handle low RNA input samples in a tube reaction (Online Methods). Next, we subjected each sample to RNA-seq, generating 1.6 million raw reads per cell, of which 40% mapped uniquely to 2,068 distinct genes (on average, using UCSC conservative gene models, mainly protein-coding). On average, we observed 5,977 and 6,760 mRNA molecules in excitatory and inhibitory neurons, respectively (Fig. 3d). This corresponds to an absolute capture efficiency of 7% per mRNA molecule per cell, as inferred by comparison with previously published single-cell data 11 . By comparison, RNA-seq of single neurons from the mouse neocortex can recover ~19,000 RNA molecules/neuron, mapped to ~5,000 distinct genes 11 , with a capture efficiency of ~20%. We attribute losses to procedural differences, leakage during cell aspiration and/or binding of aspirated RNA to glass or plastic surfaces. Pan-neuronal markers 11 Thy1 and Stmn2 were detected in 79% and 87% of the cells, respectively. To validate the quality of our RNA-seq data, we compared the results from inhibitory and excitatory neurons. As expected, I-type neurons expressed Gad1 (the gene encoding GAD67; 41 of 45 cells) and Gad2 (the gene encoding GAD65; 39 of 45 cells), as well as Cck (44 of 45 cells). Moreover, all subtypes of interneurons contained mRNAs for the vesicular inhibitory amino acid transporter (Slc32a1) 31 but lacked either vesicular glutamate transporter 1 (Slc17a7) or 2 (Slc17a6) 32,33 . Subsets of these interneurons also expressed mRNAs for neuropeptides (Vip, Npy and Crh) and Ca2+-binding proteins (Calb1; Fig. 1g, 3e), markers known to co-exist with CCK in interneurons 5,6 . None of the cells (0/45) contained parvalbumin (Pvalb) mRNA transcripts. By contrast, Exc-type pyramidal cells only expressed Slc17a7, and to a lesser extent Slc17a6, but not Slc32a1 (0 of 38 cells), along with considerably lower copy numbers of Cck mRNA (Fig. 3d). Thus, our RNA-seq data accurately reflect the major functional distinction between inhibitory and excitatory neurons in the cerebral cortex. Examining the RNA phenotypes of each of the 5 I-types amongst interneurons revealed a distinct pattern of common molecular markers (Fig. 1g, 3e). Thus, at this level of analysis, there was a one-to-one correspondence between transcriptionally- and electrophysiologically-defined cell identities. Mapping neuronal identities on single-cell RNA-seq datasets Electrophysiology is inherently limited in throughput. Consequently, the molecular classification of neurons from small and/or heterogeneous groups of cells is challenging because of the resulting low statistical power. To increase the reliability of our molecular classification, we took advantage of the single-cell dataset on somatosensory cortex we recently generated 11 , containing >3,000 single-cell transcriptomes. We reasoned that if Patch-seq data could be aligned to this much larger dataset, we would be able to assign electrophysiological properties to molecularly better-defined neuronal subclasses. Even though this approach is not mandatory for neuronal classification, the increasing availability of reference datasets for major brain regions will enhance overall classification accuracy in small-sized sample populations. We built a correlation-based classifier to assign each neuron from Patch-seq to one of the possible neuronal subtypes distinguished earlier 11 . The classifier used an iterative process of selecting relevant features (e.g., genes), ranking the candidate groups by correlation with any measured cell, and removing groups with lower correlation (Fig. 4a and Online Methods). None of the excitatory cells (0/38) were classified as interneurons, whereas a single interneuron (1/45) was classified as a L5 pyramidal cells (Fig. 4b). The anatomical position, electrophysiology classification and molecular phenotype of pyramidal cells with high certainty corroborated their layer-specific L2/3-L5 identities (Fig. 4c). All of the measured interneurons were assigned to CCK+ subclasses (Fig. 4d). None of the “I-type” interneuron subdivisions was classified into a single group molecularly. Instead, I-types 1 and 5 were aligned to one set of closely related interneurons (Int11-Int14 in ref. 11 ; Fig. 4d). Likewise, I-types 2 and 3 were assigned to Int5-Int8, forming another subset of closely-related cell types. I-type 4 showed exceptional heterogeneity with half of the cells assigned to each of these subsets. This shows that the combination of patch-clamp and RNA-seq methods benefits from an increased power of classification through the combination of real-life biophysical (and morphological) criteria and statistical predictions. Molecular candidates to determine interneuron excitability The depth of our molecular analysis allowed us to quantitatively assay the expression of channels, ion pumps and receptors in the I-type interneurons (Fig. 5a-c). Having our cells patch-clamp recorded, we could compare expression differences of any subunit detected with membrane potential changes (Fig. 5d). For example, Na+/K+ adenosine triphosphatase (ATPase) is a key electrogenic determinant of Vrest in excitable cells 34,35 . The quantitative expression of genes encoding ATPase subunits (Atp1a1-Atp1b3; Fig. 5a) closely and positively correlated with Vrest of each neuronal subtype in our dataset (for Atp1a3 see Supplementary Fig. 3a). Voltage-gated Cl− channels, a family of poorly understood ion channels 36 , are thought to modify Vrest by gating ion fluxes. Here, Clcn3 was shown to positively correlate with Vrest (Supplementary Fig. 3b), potentially implicating these channels in determining subthreshold membrane potential fluctuations. For AP frequency modulation, Kv3.1 (Kcnc1), a delayed rectifier K+-channel, is broadly assumed to confer the capacity to discharge at high frequencies upon prolonged depolarization of fast-spiking interneurons 37 . Our data support this because interneurons firing in excess of pyramidal cells expressed ~3x higher RNA count/cell for Kcnc1 (p < 0.05 for the whole I-type cell population). Similarly, coherence of the electrophysiology and RNA data was demonstrated by detecting cyclic nucleotide-regulated ion channel Hcn1. We saw no expression in “I-type 1” cells, where electrophysiology confirmed the lack of hyperpolarization-activated “sag” depolarization, a biophysical consequence of Hcn1 activity (Fig. 1b; Supplementary Table 1). Thus, our RNA-seq data will allow predictions for future neurophysiology studies interrogating specific parameters of neuronal excitability. Correlation matrix for use-dependent markers of neurons Independent of any a priori classification, Patch-seq also permitted the analysis of correlations between gene expression and electrophysiological parameters. Many genes (748 out of 5,600) showed significant correlation with one or more electrophysiology parameters. We took advantage of our quantitative datasets, and asked if rendering mRNA copy numbers of ion channels and synapse-related proteins (167 passed our criteria, Online Methods) as predictors returns meaningful association with specific biophysical parameters of single APs or AP trains. We hypothesized that any sufficiently robust correlation (filtered for correlation coefficients exceeding −0.4 and/or +0.4 11 ; Supplementary Fig 4a,b) could be valuable for future studies if it allows for distinguishing any CCK+ interneuron subtype. 24 genes (Clic4, Clip3, Cacna1g, Kcnma1, Kcnj11, Kcnc1, Apba2, Cacna1g, Cadps2, Exoc8, Gria1, Grin2b, Htr7, Kcnma1, Npy, Pak1, Pcdh8, Slc32a1, Slc6a17, Sort1, Stx4a, Syt6, Syt7, Tac2) had significant correlation with at least one parameter (Supplementary Fig. 4d-i). The main advantage of this approach is that it provides testable hypotheses by focusing on the preferential expression of genes in one but not another subset of interneurons. For example, synuclein-γ (Sncg), synaptotagmin 7 (Syt7), vesicle-associated membrane protein 4 (Vamp4) and the GABA synaptic reuptake transporter Slc6a1 were grouped together and associated with AP parameters. Although the exact functional significance of these associations remains elusive, they provide candidates for the future molecular dissection of neuronal networks under physiological or pathological states. Subtype-specific receptor repertoire in CCK interneurons At any point in time, the intrinsic excitability of a neuron is dynamically tuned by its afferent inputs. An advantage of our Patch-seq dataset is that it contains information on most (if not all) ligand-gated ionotropic channels, metabotropic (G protein-coupled) and other receptors, which determine the net network load on each sampled neuron (Fig. 5a-c). This allows for inferences be made on the specificity and heterogeneity of afferent inputs. For example, fast glutamatergic transmission relies on α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptors expressed by all neurons. Its subunits, GluR1 (Gria1)-GluR4 (Gria4), are differentially expressed in various neuronal subtypes. As such, the fast kinetics of AMPA receptor inactivation and desensitization in interneurons is usually explained by their low levels of GluR2 (Gria2) subunit 38 . Here, we sampled 210 Gria1 mRNAs in total in 45 interneurons (4.46 ± 5.3 mRNA molecule/cell) and 283 mRNAs cumulatively in 38 pyramidal cells (8.57 ± 8.08 mRNA molecule/cell; p < 0.05, mean difference of 3.61; Fig. 5b), thus quantitatively recapitulating earlier predictions by histochemistry and channel neurophysiology. Next, we examined the expression of the type 1 cannabinoid receptor (Cnr1) in our sample, which was reported to be highest in CCK+ interneurons in the brain 39,40 . Here, we find that I-type 1 cells lack appreciable Cnr1 mRNA expression (12 Cnr1 mRNA molecules in total), which contrasts with I-type 2-5 subclasses (678 Cnr1 mRNAs in total; Fig. 5c). These differences suggest another hitherto undescribed level of molecular complexity amongst cortical CCK+ interneurons. Many developmental biology studies utilize the serotonin (5-HT) receptor 3a (Htr3a) subtype to detect and classify CCK interneurons 4,41 . In our Patch-seq dataset, we observed no 5-HT3a expression in I-type 1 cells (Fig. 6a,c), low Htr3a mRNA expression in I-type 2 (Calb +/Npy + /Vip −) cells and high mRNA copy numbers in I-type 3 (Vip +/Npy −) and I-type 4 (Vip +/Npy +, Fig. 6b,c) and I-type 5 (Vip −/Npy +) interneurons. This suggests that Htr3a expression amongst CCK+ interneurons might be more restricted than initially thought, likely leading to the under-sampling of CCK interneurons by Htr3a-driven genetic tags 41 . Subsequently, we expanded our analysis to all Htr1a-Htr7 subtypes (Fig. 6c). None of the CCK+ interneurons contained Htr1a, Htr1b or Htr1d mRNAs. Notably, I-type 1 cells contained no mRNA copy for any of the 5-HT receptors. We validated our Patch-seq results by measuring the excitatory effect of 5-HT on I-type 1 (n = 4) vs. I-type 4 (n = 4) interneurons (Fig. 6a,b) in whole-cell current-clamp experiments. While clamping the cells to 0 pA, 5 μM (data not shown) and 25 μM 5-HT (Fig. 6d,e upper trace) was bath-applied. Both subthreshold depolarization with 5 μM and suprathreshold depolarization with 25 μM 5-HT 42-44 was observed in both I-type 1 and 4 interneurons. To dissect the origin of the depolarizing 5-HT effect and rule out indirect effects imposed by the neuronal network the recorded cell was embedded into, 10 μM 5-HT was microejected (”puffed”) onto interneuron somata. We controlled equal 5-HT load by coapplying a fluorescent tracer (Fig. 6d,e middle). Direct 5-HT application onto I-type 1 interneurons did not cause depolarization (Fig. 6d, bottom). By contrast, I-type 4 interneurons became readily depolarized (Fig. 6e, bottom). Our pharmacological rsults thus establish causality, considering that RNA-seq predicts 5-HT insensitivity in I-type 1 interneurons due to the lack of any metabortopic or ionotropic 5-HT receptors. Moreover, the morphological reconstruction of the sampled interneurons (Fig 6a,b) is explanatory towards the depolarization of I-type 1 cells upon bath application of 5-HT: I-type 1 interneurons were dye-coupled to neighboring cortical non-pyramidal neurons, and even pyramidal cells. This allowed for interneurons to synchronize their membrane potentials through gap junctions 45 , assigning I-type 1 cells as ”passive followers” of network depolarization at suprathreshold 5-HT load. Thus, our Patch-seq data primed us to identify a CCK+ interneuron with a 5-HT-driven effector behavior in its L1 microcircuit even if itself lacked 5-HT receptors. Overall, we suggest that our Patch-seq approach will help future systems (neuro-)biology investigations to rationalize diverse functional outcomes by providing a platform of identifying marks for specific modalities amongst neurons. Discussion The brain undoubtedly exhibits the highest level of cellular heterogeneity 10 , and contains an large variety of neurons that differ in their morphology, connectivity, biophysical parameters and molecular phenotypes 5,6,46 . The taxonomy for neurons dates back to the first pioneers of neuroanatomy (e.g., Cajal and Golgi), who exclusively used morphological features, such as the size and topography of axonal and dendritic arbors, for classification and is now based on a wide array of neurophysiology, advanced histochemistry and RNA analyses 6 . Nevertheless, reliance on known candidate marks continued to dominate, and delayed the inception of unbiased classification. Moreover, the limited number of markers that could be probed at any given time (~20 for single-cell PCR and histochemistry) 7 together with the often mutually-exclusive experimental conditions that neurophysiology and single-cell molecular biology tools require limited detailed fingerprinting of cellular components in the brain. We combined mouse genetics and patch-clamp electrophysiology to successfully target a distinct cohort of interneurons 5,47 in order to overcome existing limitations of classical function-structure analyses. In the CCK+ cell population studied, we expected considerable neuronal diversity, which allowed us to improve and optimize Patch-seq to its present accuracy, even when low copy numbers of mRNA molecules were present. We first classified our CCK+ interneuron sample from L1/L2 (but not from deep cortical layers 11 or hippocampal subfields 5,47 ) into the 5 commonly identified subtypes. We further showed that Patch-seq can identify a select number of molecular determinants that can be used to further subdivide CCK+ interneuron subclasses. Each of these examples is significant because they reconcile previously reported sets of data on receptors (5Ht3a, Cnr1) and channels (Hcn1). We also note that Patch-seq can discover sets of cellular markers independent of other classification systems or a priori knowledge of the cell-type of interest . We also show that Patch-seq can help to form hypotheses e.g. about neurotransmitter-receptor relationships. We expect that similar approaches can be applied to essentially any neuronal subtype, and will help to avoid false-negative data (due to under-sampling of neuronal contingents) in cases when the resolution of available histochemical tools is limited or if suitable reagents are not available. As such, the lack of Rgs12, Htr3a and Cnr1 mRNA transcripts in I-type 1 interneurons raises the possibility that this abundant subclass might have been systematically missed in prior genetic reporter analyses, curtailing the analysis of their contribution to fundamental cortical network events. The efficiency of mRNA capture in Patch-seq is lower than that in single-cell RNA-seq on dissociated tissues. However, it is still sufficient to efficiently sample even “low” expressed genes because of its extremely low rate of false-positive identification 11,12 . Thus, mRNA copy numbers even at the range of 1-5 molecules return meaningful associations. Moreover, the combination of Patch-seq with transgenic mouse technologies might allow the future exploitation of external (“spike-in”) reference standards (e.g. transgene products), thus facilitating positive cell identification. These methods together with the progressive decoding of regional heterogeneity in the nervous system through large-scale RNA-seq databases 11-13 can increase the stringency of neuronal classification. Such “reference atlases”, once available, will allow for precise hierarchical landscapes be built even when cell numbers from patch-clamp electrophysiology experiments are limited. However, Patch-seq can stand alone and give much more complete and accurate information about gene expression (~2,000 genes per cell) in selectively probed cell contingents, compared to previous methods (e.g. qPCR for 10-20 genes/cell) 6,7 . Patch-seq samples somatic material upon aspiration. For neurons, dendrites and axons occupy large spaces and their intracellular volume is considered significant. Therefore, one may argue that Patch-seq misses many mRNAs that are preferentially targeted to distant domains of axons or dendrites. Although some mRNA is actively transported into neurites 48 , this does not mean that they are absent from the soma. On the contrary, most (if not all) mRNA species are more abundant in the soma than in neurites, and there is not a single known case of an mRNA that is localized exclusively outside the soma. This is true even for mRNAs thought to be actively transported into neurites, such as CamKIIα and spinophilin 49 . It is also known that the axon and dendrite contain much less total mRNA than the soma: axons are thought to contain about 1,000 – 4,500 mRNAs, whereas dendrites contain >2,500 mRNAs 49 . These numbers are at least an order of magnitude lower than the soma (which contains >100,000 mRNAs). As a result, sequencing the soma content can be expected to give a representative view of mRNA expressed by a neuron, although without information on which mRNAs are more or less efficiently transported into the neurites. Another technical element that needs to be tightly controlled is the length of electrophysiology recordings because electrical stimuli might alter the transcriptome. Here, we used 20-25 min protocols, which are compatible with the life-time of mRNA being on the order of many hours (median: 9 hours) with none being known as shorter than 1 hour 50 . The quickest transcriptional response known in any setting is the induction of immediate early genes (e.g., c-Fos, Jun, Egr2), and can be detectable after 30 min but peaks at 3 hours. However, removing immediate early genes did not affect cell-type identification. Thus, the impact of patch-clamp recordings that occur on a timescale less than one hour can be expected to have minimal impact on the RNA transcriptome. In conclusion, Patch-seq can be expected to facilitate the characterization of transcriptome-wide changes in many experimental settings, thus contributing to a better understanding of fundamental physiological and pathological processes. Online methods Animals and husbandry We generated a dual reporter mouse line (Supplementary Fig. 1a) by crossing parental lines that either expressed Discosoma red fluorescent protein (DsRed) under regulatory elements of the cholecystokinin (CCK) promoter on a bacterial artificial chromosome (BAC/CCK-DsRed) 24 or green fluorescent protein (GFP) knocked into the glutamate decarboxylase 67 gene (GAD67+/gfp) 23 . The resulting CCKBAC/DsRed::GAD67+/gfp line appeared anatomically normal, particularly without changes to brain size or deformities to its fine structures, including normal cell proliferation, migration (data not shown) and laminar distribution (Supplementary Fig. 1a). Animals were group-housed under 12:12 h/h light/dark cycle with ad libitum access to water and food. Animals of both sexes were used for neurophysiology experiments during postnatal days 17-23. Experiments on live animals conformed to the 86/609/EEC directive and were approved by the regional authorities on animal ethics (Stockholm Norra Djuretiska Nämnd; N512/12; Tierversuchsgesetz 2012, BGBI, Nr. 114/2012). Preparation of brain slices, correlated differential-interference contrast and epifluorescence microscopy, superfusion All experiments on interneurons were performed in L1/2 of the primary somatosensory cortex (S1). Coronal slices (300-μm thickness) were prepared on a VT1200S vibratome (Leica, Germany) in ice-cold artificial cerebrospinal fluid containing (in mM): 90 NaCl, 2.5 KCl, 1.25 Na2HPO4, 0.5 CaCl2, 8 MgSO4, 26 NaHCO3, 20 D-glucose, 10 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES), 3 Na-pyruvate, 5 Na-ascorbate (pH 7.4). Brain slices were then incubated at 22-24 °C for 60 min in a “recording” solution containing (in mM): 124 NaCl, 2.5 KCl, 1.25 Na2HPO4, 2 CaCl2, 2 MgSO4, 26 NaHCO3, 10 D-glucose (pH 7.4). All constituents were from Sigma-Aldrich. Both solutions were aerated with carbogen (5% CO2/95% O2). Temperature was set to 33 °C (TC-10, Npi, Germany) in the recording chamber. Brain slices were superfused with “recording” solution at a rate of 4-6 ml/min. Neurons were visualized by differential interference contrast microscopy on an Olympus BX51WI microscope. Next, the co-existence of DsRed and GFP in interneurons was confirmed by epifluorescence microscopy using a mercury arc light source (USH-1030L, USHIO) and appropriate combinations of band- and long-pass excitation/emission filters (for GFP: U-MWIG3 [excitation: 460-495 nm; emission: 510-550 nm]; for DsRed: U-MWIBA3 [excitation: 530-550 nm; emission: >575 nm]). Serotonin (5-HT; Tocris) was directly dissolved in the recording solution at final concentrations of 5 μM, 10 μM and 25 μM, and superfused at flow rates as above. Focal 5-HT ejection was performed using a microinjector (PDES-02TX-LA, Npi, Germany) after filling borosilicate glass capillaries (Hilgenberg, 3-4 MΩ) with 10 μM 5-HT. Pressure pulses of 500 mbar for 0.5 s were used with 30 s intervals. Pharmacological probing of the interneurons was carried out at 33°C. Patch-clamp electrophysiology Recordings were carried out with borosilicate glass electrodes (Hilgenberg, Germany) of 3-4 MΩ pulled on a P-1000 instrument (Sutter, USA). Electrodes were filled with an intracellular solution containing (in mM): 130 K-gluconate, 6 NaCl, 4 ATP-Na2, 0.35 GTP-Na2, 8 phosphocreatine-Na2, 10 HEPES, 0.5 ethyleneglycolbis(2-aminoethylether)-N,N,N',N'-tetraacetate (EGTA) and 0.5 mg/ml biocytin (pH 7.2 set with KOH). Whole-cell patch-clamp recordings were made on an EPC-10 triple amplifier (HEKA, Germany). Current clamp recordings were corrected for −9.99 ± 0.38 mV liquid junction potential between the intracellular and recording solutions, as measured against a 3M KCl-electrode. Resting membrane potential (Vrest, expressed as mVs) was measured in current-clamp mode at 0 pA current. Input resistance (Rm, expressed as MΩ) was calculated using linear regression established between electrotonic voltage responses (±15 mV from Vrest) and 500-ms current steps of increasing amplitude (10 pA increments). Membrane time constant (τ, ms) was averaged from 20 successive electrotonic voltage responses to hyperpolarizing (−40 pA) current steps and fitted with a single exponential. The “sag” depolarization, indicating the activation of hyperpolarization-activated non-selective cationic current (I h ) was calculated as the ratio of the peak negative voltage and steady-state negative voltage at a voltage response to hyperpolarizing current injections resulting in a steady-state voltage of −100 mV. AP threshold (APthr, mV) was defined as the voltage point where the upstroke’s slope trajectory first reached 10 mV/ms. AP amplitude was defined as the difference in membrane potential between APthr and AP at peak. AP rise time (ms) was the time from the APthr to the AP’s peak. AP duration (ms) was the time interval from APthr to the same voltage value during repolarization. The amplitude of afterhyperpolarization (AHP, mV) was defined as the difference between APthr and the most negative membrane potential attained during the AHP. AP decay time was calculated as the time from the AP’s peak to the peak of the AHP. Afterdepolarization (ADP) amplitude (mV) was defined as the difference between AHP peak and the most positive voltage value between AHP peak and the fast repolarization peak. In case of a lack of AHP, ADP was recognized as the most positive membrane potential between steady-state voltage and the peak of fast repolarization. Maximum AP up- and AP down-stroke were determined as the maximum and minimum of the geometrical differential of the AP (mV/ms), respectively. Maximum up- and down-stroke times were the times from APthr to reach maximum AP up- and down-stroke, respectively. These parameters were measured for i) the first AP elicited by 500-ms rheobasic current step; ii) the first AP evoked along a 1-s current ramp of 0-150 pA and iii) the first three APs generated on a 750-ms trace as the voltage deflection of double of the rheobasic current was injected. Adaptation ratio was calculated as the ratio of the last inter-spike interval relative to the first five inter-spike intervals. Firing frequency (Hz) was determined at saturating current injections producing spike trains. AP amplitude accommodation was determined as a ratio of the average of the last three AP peaks relative to the first five AP peaks (difference defined in mVs). All parameters were measured (Supplementary Table 1) by applying manual procedures custom-written in Matlab (MathWorks, USA). Tissue clearing and light microscopy Brain slices containing biocytin-filled neurons were post-fixed in 4% paraformaldehyde in phosphate-buffer (PB, 0.1M, pH 7.8) at 4 °C overnight. Slices were repeatedly washed in PB and cleared using “CUBIC reagent 1” (25 wt% urea, 25 wt% N,N,N’,N’-tetrakis(2-hydroxypropyl) ethylenediamine and 15 wt% polyethylene glycol mono-p-isooctylphenyl ether/Triton X-100) for 2 days 25,26 . After repeated washes in PB, biocytin localization was visualized using streptavidin-conjugated horseradish peroxidase (Vectastain Elite) using 3.3′-diaminobenzidine (DAB) as chromogen and H2O2 (0.05%) as substrate (in Tris-HCl, pH8.0). Slices were then re-washed in PB and submerged in “CUBIC reagent 2” (50 wt% sucrose, 25 wt% urea, 10 wt% 2,20,20’-nitrilotriethanol and 0.1% v/v% Triton X-100) for further clearing. Post-hoc neuroanatomy and neuronal reconstructions were performed using 40x and 60x objectives on a light microscope (Olympus BX51). Optionally, cleared tissues were imaged on a laser-scanning microscope (LSM780 and ZEN2013 software, Zeiss). Three-dimensional filaments of DAB-stained cells from 300 μm slices were reconstructed in Neurolucida (cx9000, Mbf Bioscience). Cell harvesting for sequencing A the end of each patch-clamp protocol, the micropipette was clamped to a holding potential (Vhold) of −5 mV. Prior to the harvesting procedure, a continuous series of depolarizing rectangular voltage pulses (5 ms at 5 ms intervals) were applied for 6-7 min with amplitudes of 25 mV from Vhold. The entire soma of each recorded neuron was aspirated into the micropipette slowly (~1-2 min) by applying mild negative pressure (−50 mPa). This procedure allowed us to retain a tight seal and to minimize RNA loss by keeping the RNA molecules in the pipette solution. When we broke contact, the recording pipette was pulled out from the recording chamber and then carefully rotated over an expelling 0.2 μl tube, where its content (0.8-0.9 μl) was ejected onto a 0.6-μl drop of lysis buffer pre-placed onto the side of a 0.2 ml tight-lock tube (TubeOne). The resultant sample (1.5 μl) was span down (20 s) at 24 °C to the bottom of the conical tube, stored at −80 °C and later subjected to in-tube reverse transcription (RT). Lysis, cDNA synthesis and library preparation Cell aspirates were dispensed into ~0.5 μl lysis mix consisting of 0.15% Triton X-100, 1 U/μl TaKaRa RNase inhibitor, 4 μM reverse transcription primer C1-P1-T31 5′-Bio-AATGATACGGCGACCACCGATCG TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT-3′, 3.5 mM dNTP and 17.5 mM DTT. Samples were collected and stored at −80 °C until batch processing. Before reverse transcription, samples were thawed and lysed at 72 °C for 3 min, then cooled to 4 °C. Immediately following the lysis step, 2 μl RT mix (1x SuperScript II First-Strand Buffer; Life Technologies) supplemented with 10.6 mM MgCl2, 3.6 μM template-switching oligo C1-P1-RNA-TSO 5′-Bio-AAUGAUACGGCGACCACCGAUNNNNNNGGG-3′, 1.5 U/μl TaKaRa RNase inhibitor (Clontech), 1.45 M betaine and 21 U/μl Superscript II reverse transcriptase (Life Technologies)) were added and incubated at 42 °C for 90 min followed by 72 °C for 10 min. Following reverse transcription, 8 μl PCR mix (1x KAPA HiFi 2x ready mix and 240 nM C1-P1-PCR2 5′-Bio-GAATGATACGGCGACCACCGAT-3′) were added and PCR-amplified using thermal cycling as follows: 95 °C for 3 min (5 cycles), 98 °C 20 s, 62 °C 4 min, 72 °C 6 min, (9 cycles) 98 °C 20 s, 68 °C 30 s, 72 °C 6 min, (7 cycles) 98 °C, 30 s, 68 °C 30 s, 72 °C 7 min. Subsequently, PCR samples were cleaned using AMPure-XP beads (1:1 ratio; Beckman Coulter) and quantified by qubit (Life Technologies) on an Agilent bioanalyzer. Library preparation was done using tagmentation as described 29 . Illumina sequencing Libraries were sequenced on an Illumina Hiseq2000 instrument using C1-P1-PCR2 as read1 (50 nt) primer and C1-TN5-U PHO-CTGTCTCTTATACACATCTGACGC as index read (8 nt) primer. Data analysis Read processing and molecule counts were performed as reported recently 11 . We only analyzed cells with >1,500 mRNA molecules/cell (excluding mitochondrial, repeat and rRNA) and if a complete catalogue of patch-clamp read-outs was available. Alignment of interneurons and pyramidal cells on a cortical template We used our recently described cortical dataset 11 to resolve each of the interneurons and pyramidal cells into one of template groups (from >3,000 dissociated cells). First, we narrowed down our search to one of the layer-specific pyramidal cell cohorts (S1PyrL1-L6) or interneuron groups (Int1-16), 22 groups in total. Because of a significant difference in the number of mRNA molecules detected, we designed our classifier on correlation measures rather than Euclidian distance. As feature selection is an important parameter for classification (e.g., which genes are “in use”), the classifier continuously updated the features for the groups compared. First, the median expression for every group (in the cortex dataset) was calculated. Because the standard deviation of the genes’ median expression (along groups) did not depend strongly on their mean expression, we selected genes (features) using a fixed threshold of standard deviation being >1.5 over the groups compared. The process was as follows: i) selection of a Patch-seq cell for classification; ii) setting the list of candidate groups to contain all 47 neuronal groups of the somatosensory cortex; iii) selection of genes for the current iteration, std(gene median expression) >1.5; iv) calculation of the correlation between a Patch-seq cell and all candidate groups; v) ranking neuronal groups by their correlation to the Patch-seq cell (highest to lowest); vi) removal of those candidate groups that were <50% with their correlation; vii) if the list of candidate groups contained more than one group we looped back to step iii. If only a single candidate group was left, the Patch-seq cell was assigned to that particular group. Correlation of gene expression and electrophysiology parameters This analysis was aimed to identify mRNAs that might be predictors for electrophysiology parameters. We focused on mRNAs coding for ion channels, receptors and synaptic transmission-related proteins as sets of genes whose biological interpretation might be tested experimentally. We tested the correlation between all genes (5,600) that passed our baseline criteria (that is, more than 5 cells with non-zero expression, ~12% of data points) against all electrophysiology parameters (110). Along with the correlation, p-values for the null hypothesis of independent variables were calculated (permutation p-value shows similar results). For each gene separately, we used a false recovery rate (FDR) of 10% to declare a significant correlation. This was performed separately for every gene because otherwise the p-values strongly correlated and the assumptions for FDR were violated. Next, we focused only on ion channels, receptors and synapse-related genes (167 of these passed our above criteria). Correlation coefficients shown in Supplementary Fig. 4 are for genes that exhibited correlation (or anti-correlation) greater than 0.4. Supplementary Material Supplementary Figures 1-4 Supplementary Table 1