Single-cell RNAseq reveals cell adhesion molecule profiles in electrophysiologically defined neurons

As a starting point, we examined the transcriptome of the entire hippocampus (Fig. 1A) at five developmental stages (in triplicate at postnatal days P0, P7, P14, P21, and P28). In these samples, the total number and distribution of expressed genes were similar (Fig. 1 B and C, and Fig. S1). To specifically analyze cell adhesion molecules, we created a comprehensive list of candidate genes involved in, or related to, transsynaptic cell adhesion (collectively referred to as CAMs, for candidate cell adhesion-related molecules). For this, we first considered all molecules with single transmembrane domains (a general, although not unique property of membrane surface adhesion molecules) and narrowed this list down to 406 genes based on preexisting data (Fig. S1C). This curated list of candidate cell adhesion molecules includes proteins implicated in cell–cell signaling but excludes, for example, intracellular signal transducers linked to cell adhesion signaling. We found these genes to be broadly represented in all replicates and developmental stages (Fig. 1D), highlighting diversity of CAMs in the hippocampal circuit and throughout development. Expression profiles (Fig. 1E) for each developmental stage included genes that were consistently highly expressed at all developmental stages [for example, amyloid precursor protein (App), contactins (Cntn), neurexins (Nrxn), and protein-tyrosine phosphatases (Ptpr); Fig. S1D], as well as genes with large changes during development (for example, Ncam1 and Sdc3). Normalization of individual gene expression levels also revealed an average change of 60% between P0 and P28 (Fig. 1 F and G).

To analyze gene expression in specific cell types in a defined circuit, we developed a method—which was also recently introduced by others (16, 17)—with which we first characterized neurons electrophysiologically, and subsequently analyzed their transcriptome by aspiration of cytosol followed by single-cell RNAseq. Using this approach, we examined FS and RS interneurons that directly inhibit CA1-PYR cells, as well as the target CA1-PYR cells for these interneurons. We patched FS and regular-firing interneurons located within, or in close proximity to, the pyramidal cell layer in wild-type mice, and then used paired recordings by simultaneously patching PYR cells to test whether presynaptic action potentials (APs) evoked inhibitory postsynaptic responses, characterizing these cells as FS or RS inhibitory interneurons (Fig. 2A). Note that such distinction did not necessarily identify individual cell types. For instance, FS cells could include PV+ basket, bistratified, and axoaxonic cells, whereas RS cells could include CCK+ basket or bistratified cells as well. At the end of the recordings, we collected the neuronal cytosol via pipette aspiration for RNAseq (Fig. 2A). Upon alignment of sequencing reads and assignment of gene counts for each gene, we applied stringent quality control metrics to each cell (Fig. S2 A and B). RNAseq data obtained in this manner from 18 RS-INT, 9 FS-INT, and 14 CA1-PYR cells passed quality control.

To gain insight into the molecular identity of recorded cells, we examined expression of genes that had been previously associated with RS-INT, FS-INT, or PYR cells (Fig. 2B). As expected, Gad1 and Gad2 were highly expressed in interneurons but not in PYR cells, consistent with their respective neurotransmitters. Although CCK peptide is a unique marker for RS CCK cells, we found that the CCK gene was nonspecifically expressed in all three cell types. Similarly, although cannabinoid type-1 receptor (Cnr1) is most highly expressed in RS CCK cells (for specific examples, see refs. 10, 18, and 19; for review, see ref. 7), it was also produced in FS-PV and PYR cells. In contrast, Htr3a and PV (parvalbumin) were more specific to RS-INT and FS-INT cells, respectively. Although the latter is a unique marker for FS subtypes at the protein level, the corresponding mRNA was detected in two RS-INT and three PYR cells as well. To identify PYR cells, we examined Camk2a and Neurod6 expression: although the former’s expression was rather nonspecific, Neurod6 was a reliable predictor of PYR cell identity. We also found that Grik3 (kainate receptor GluR7) uniquely, but not exclusively, identified FS-INT cells.

Examining overall gene expression, we found that interneurons consistently expressed almost twice as many genes as CA1-PYR cells (Fig. 2C; see Fig. S2F for controls), and that detection of individual genes was also more consistent in RS- and FS-INT cells than in CA1-PYR cells (Fig. 2D). A potential explanation for the latter observation involves spatial gene expression gradients in the hippocampus that likely underlie heterogeneity of CA1 pyramidal neurons (20).

To link single-cell transcriptional profiles to functional properties, we examined expression of 140 voltage-gated ion channels (see Materials and Methods for gene set assembly; Fig. 3A, Fig. S3A, and Dataset S1). Because RS-, FS-, and CA1-PYR cells have distinct electrophysiological properties (Fig. 3A), we asked whether ion channel expression could account for their physiological differences. On average, we detected more ion channels and more consistent gene expression in interneurons than in PYR cells (Fig. 3 B and C), following total transcript levels. As a pivotal example of cell type-specific gene expression, we found enrichment of Na⁺ and K⁺ channels in FS cells (Fig. 3 D and E, Upper).

Next, we quantified electrophysiological parameters (n = 11; Fig. 3 A and E, Left) because their correlation with gene expression patterns could reflect characteristic, cell type-specific properties in these cells. We used unsupervised clustering and plotted genes/electrophysiology parameters with at least one correlation value of greater than 0.35 or less than −0.35 (arbitrary threshold to eliminate no or minimal correlations; Fig. 3E). Here, we found that AP firing rate, threshold, symmetry, and afterhyperpolarization, parameters that typically exhibit high values in FS cells, correlated with high expression of Na⁺ and K⁺ channels (Fig. 3E), including Kcnc1 (Kv3.1) and Kcnc2 (Kv3.2). Furthermore, we found that expression of these genes inversely correlated with AP peak-to-trough and width times (i.e., larger gene expression values correspond to lower electrophysiological property values, and vice versa), which as properties also correspond to the fast signaling characteristics of FS-INT cells. These results are not unexpected as expression of Kcnc1 and Kcnc2 is a hallmark of FS PV interneurons (refs. 21 and 22; for review, see ref. 23), validating the single-cell RNAseq analysis.

Surprisingly, unbiased correlation analysis can also deliver seemingly contradictory results. For example, the observation of high Hcn1 expression but a low hyperpolarization sag in FS-INT is puzzling because a primary readout for Hcn-mediated I_h currents is the sag amplitude. However, the correlation between hyperpolarization sag and I_h current is not absolute, nor the involvement of Hcn1 in modulating active membrane properties (24). In addition, we found significantly higher expression of Kcnc3, Scn1a, and Trpc6 (Fig. 3E) in FS cells (described for PV cells in refs. 25–27, respectively), whereas specific expression of Kcng4 (Kv6.3, Kv6.4) identifies a thus far uncharacterized component of the ion channel repertoire of FS cells. For further clues about molecular architecture of these cell types, see expression of ligand-gated ion channels in Fig. S3 C–E.

Together, these combined analyses extend previous single-cell studies (including refs. 22 and 26–29) that used probe-based gene expression analysis, and enable comprehensive and unbiased RNAseq-based functional inferences.

Our analyses in whole-tissue preparations (Fig. 1) consistently detected transcripts encoding 383.2 ± 1.1 (mean ± SEM) different CAMs in the hippocampus, which comprises a large variety of neuronal and nonneuronal cell types. In examining expression of CAMs at the single-cell level, we focused on how many CAMs were expressed in a neuron, and whether these were expressed with cell type specificity. In single RS-INT, FS-INT, and CA1-PYR cells, we detected N_RS-INT = 121.6 ± 11.2, N_FS-INT = 128.7 ± 12.7, and N_CA1-PYR = 82.7 ± 4.8 different CAM transcripts, respectively [P_{RS-INT vs. FS-INT} = 0.73, P_{RS-INT vs. CA1-PYR} = 0.02, P_{FS-INT vs. CA1-PYR} = 0.006; pairwise Mann–Whitney (MW) test; Fig. 4 A and B and Dataset S1]. Examining these genes revealed that expression of CAMs partially overlapped between different cell types. There were N_{RS-INT vs. FS-INT} = 48, N_{RS-INT vs. CA1-PYR} = 58, and N_{FS-INT vs. CA1-PYR} = 62 differentially expressed CAMs between cell types (pairwise MW test, P < 0.05 in all cases). Thirty CAMs were expressed higher (n = 26; MW test, P < 0.05 in all cases) or lower (n = 4; Epha4, Mdga1, Pcdhgc5, Sema3e; MW test, P < 0.05 in all cases) in both interneuron types compared with PYRs. We found that genes that were highly expressed in tissue samples were also highly expressed in the RS-INT, FS-INT, and CA1-PYR cells (note that for sequencing libraries, tissue samples were diluted and processed using same conditions and reagents as for single cells, albeit with larger amounts of starting mRNAs). Such genes included App, Chl1, Cntn1, Nptn, Nrxn1, and Ptprs that were consistently detected in all three cell types (Fig. 4C, “All cell types”).

The 26 CAMs that were expressed higher in interneurons than in CA1-PYR cells included Clstn3, Cntnap4, and Lrrc4, which have been implicated in specification of inhibitory synapses (30–32), as well as Nrxn3, a presumed presynaptic organizer of synapse function (13, 14, 33). Note that, although expression of Nrxn3 was not specific for interneurons, it was consistently enriched in interneurons, which we also considered here as a cell type-specific feature. Other examples for cell type specificity included Cbln2 and Ephb2 for RS-INT cells, Nphx1 for FS-INT cells, and Epha4 for CA1-PYR cells. Although some of these molecules have been directly implicated in transsynaptic signaling (e.g., Nrxn1, Nrxn3, Ptprs) or in synapse specialization (see examples above), others have not been studied in detail. Nevertheless, these data suggest cell type-specific expression of CAMs can potentially explain neuronal connectivity at the molecular level. Such cell type-specific features were further highlighted by principal-component axes (PCA) analysis, which revealed two main components, including overall expression as well as differential expression in interneurons vs. PYR cells (Fig. S4 A–C).

Beyond cell type-specific control of gene expression levels, diversity in CAM signaling could have been generated by other cellular mechanisms, especially by alternative splicing. Some of the most notable examples for alternative splicing include CAMs, such as Dscam in Drosophila and neurexins in mammals, which can express thousands of alternatively spliced isoforms (34–37). Therefore, we examined alternative exon use of CAMs at the single-cell level. We analyzed genes with reliable end-to-end exon junction coverage, which was applicable to n = 139 CAM genes. Out of these, we identified a single isoform in n = 67 and multiple isoforms in n = 72 genes (Fig. S4 D and E). Because of their potential cell type specificity, presence of different isoforms in the latter group would be especially interesting. However, currently available single-cell RNAseq platforms—including ours—require library fragmentation. As a result, precise exon use of full-length mRNAs can be inferred only indirectly, with sample sizes typically larger than those of single cells, hindering further analyses of cell type-specific splice variants. Nevertheless, we examined whether any splice variants of the highly expressed ubiquitous CAMs had the potential to be expressed with cell type specificity. We found that some of these molecules had only single isoforms that may or may not involve alternative splicing (for example, in Ptprs, we consistently detected exon skipping, whereas in Ptprn and Cadm, we only detected canonical isoforms), whereas others had multiple isoforms (for example, in Nptn, Nrxn1, and Cd47; Fig. S4F). Together, these exon junction observations suggested further diversification of cell adhesion signaling in single cells, which—at least regarding neurexins—was also suggested in our previous findings (9).

One presumed role of CAMs is to convey transsynaptic information and thereby establish a functional match between presynaptic and postsynaptic specializations. However, neither the precise molecular targets nor the underlying mechanisms have been described for most CAMs. Thus, we analyzed expression of two broadly defined groups of genes that are likely linked to cell adhesion function, namely genes related to vesicle exocytosis and to RhoGAP/RhoGEF signaling (see detailed list in Fig. S5A and Dataset S1). Exocytosis-related molecules (38) were expressed broadly in all single cells and at all developmental stages (Fig. 5 A and B, and Fig. S5B). In single cells, we detected both ubiquitous and cell type-specific features, where the former included Snap25, Synaptobrevin-2 (Vamp2), Snap47, and Syt11 (for which no function has been uncovered yet, despite its high abundance), whereas Cplx1, Stx1a, Syt13, and Synaptobrevin-1 (Vamp1) were expressed higher in interneurons (MW test, P < 0.05, for all genes; expression of Synaptobrevin-1 also appeared to be exclusive, as it was not detected in any PYRs; Fig. 5 A and B, and Dataset S1). In addition, we found exclusive expression of Sncg and Syt6 genes in RS-INT cells, with unknown functional consequences.

Next, we analyzed the expression of RhoGAP- and RhoGEF-domain–containing and related proteins. These are membrane-associated proteins that are thought to transform extracellular receptor-mediated signals into an intracellular response, prominently the organization of the actin cytoskeleton. Among others, the Rho/rac/CDC42 signaling machinery is thought to be involved in activity-dependent changes in postsynaptic spines (39–41). Therefore, it is plausible that RhoGAPs and RhoGEFs act as downstream effectors for cell adhesion signals. We observed diverse expression patterns for RhoGAPs and RhoGEFs between RS-INT, FS-INT, and CA1-PYR cells (Fig. 5 C and D, and Fig. S5 C and D, and Dataset S1). Particularly striking was the selectively high expression of chimaerin-1 (Chn1) in CA1-PYR cells, consistent with a possible spine-specific function (42).

Thus, far, our data offer a comprehensive view of gene expression and identified multiple cell type-specific differences. Next, we aimed to evaluate correlations in gene expression, because our data could potentially reveal genes that were developmentally coregulated and coexpressed at the single-cell level, and because such genes may contribute to core synaptic features. For this analysis, we included data on genes related to CAMs, exocytosis, and RhoGAPs/RhoGEFs, collected from five developmental stages (hippocampal tissue) and three cell types (single hippocampal cells monitored at ∼P21). To avoid inferences from inconsistently detected molecules, we excluded genes that were detected in less than five tissue or less than five single-cell samples, narrowing our analyses to 632 and 428 genes, respectively.

First, we independently examined gene correlations in tissue and single-cell data, and used unsupervised clustering to consolidate these results (Fig. 6 A and B), which reflected developmental and cell type-specific characteristics, respectively (Fig. S6A). We reasoned that functionally relevant correlations may be present in both sets with high correlation coefficients, and therefore selected gene pair correlations that were in the tail of distributions (>92nd percentile; note that this threshold was chosen arbitrarily but that consequences of this analysis were thoroughly tested at multiple threshold values; Fig. 6E) both for tissue and single-cell data (Fig. 6C). This step narrowed our focus to 269 genes and 692 correlations. Genes were outnumbered by correlations because some genes correlated with multiple others (i.e., they represented correlation “hubs”). We used graph analysis to test such possibility. In such graphs, each node represents a gene and each vertex represents correlation between two genes, if any, with a length inversely proportional to the respective coefficient; genes with more correlations simply had more vertices connecting them.

The structure of the resulting graph was surprising because two independent subgraphs emerged that exhibited dense correlations within, but no correlations between them (Fig. 6D; note the complete lack of interconnecting vertices; see Fig. S6C for complete gene listings in this graph). A biologically relevant explanation for these subgraphs was not immediately obvious, prompting us to perform additional analyses. First, we examined whether the two subgraphs were defined by nonoverlapping gene families (for example, they only include exocytosis or CAM genes). However, this was not the case, as both subgraphs included genes related to CAMs, exocytosis, and RhoGAP/RhoGEF (Fig. S6C). Second, we examined robustness of their independence by quantifying vertices between the two subgraphs while relaxing criteria for gene inclusion (from >90th to 0 percentile, i.e., all inclusion, in steps of 10 percentile). These analyses consistently revealed more intranetwork than internetwork correlations, which reversed when lowering inclusion threshold to 20 percentiles (0.2 correlation coefficient threshold in Fig. 6E, Left), suggesting robust independence (note that with “all inclusion,” 0 percentile, the graph is perfectly connected). Third, we examined randomness in graph connectivity. In a random graph, each gene would have a similar number of correlations, or vertices in the graph. However, both subgraphs displayed a nonrandom, “scale-free” nature where some genes were more interconnected than others, following a power law distribution and thus suggesting a biological origin (Fig. 6E, Right) (43). Fourth, we analyzed gene expression levels in each subgraph. Here, we found that developmental gene expression profiles were very similar within each, but not when the two sets were compared with each other. Specifically, subgraph 1 included genes that were highly expressed during early development, whereas subgraph 2 contained genes that exhibited later expression onsets, with peak expressions occurring >2 wk after birth (Fig. 6F, Left). In agreement with this finding, late expressing genes (subgraph 2) displayed higher expression than early expressing genes (subgraph 1) in single cells of all three types (Fig. 6F, Right), which were collected at ∼P21. Together, these analyses suggested developmental coregulation as a primary determinant for separation of the two graphs.

Next, we sought to further examine graph structures using added knowledge about identity of single-cell samples. We reasoned that robust gene expression correlations should be consistently present if cell types were analyzed independently. Therefore, we repeated analyses such that Fig. 6B only included RS-INT, or FS-INT, or CA1-PYR cells separately, or RS-INT and FS-INT cells combined, because they both represented GABAergic interneurons. We found that each of these analyses confirmed the independence of subgraphs 1 and 2 (Fig. S6B), and that the gene representation in each cell type-specific dataset was at least 80% identical to that in Fig. 6D (Fig. S6B), suggesting that coexpression of these genes occurs at the single-cell level. Moreover, we identified core parts of these subgraphs (i.e., parts that can be consistently detected in each cell type) by taking the intersection of the cell type-specific analyses (Fig. 6G). Because of their robust presence in each cell type, these correlations conceivably represented ubiquitous features, at different developmental stages of synapse maturation. This hypothesis was supported by the existence of correlations between structurally relevant RhoGEF and RhoGAP molecules in early development, when synapses are formed (Fig. 6G, Left, and Fig. S6D, Upper). Conversely, emerging correlations between cell adhesion and exocytosis in the late-gene network may reflect synapse specialization (Fig. 6G, Right, and Fig. S6D, Lower). Although our analyses did not reveal a mechanistic explanation for the observed correlations, they identified synaptic genes that were both developmentally coregulated and coexpressed at the single-cell level.

Our findings thus far revealed ubiquitous and cell type-specific features of CAM and signaling molecule expression as well as the contextual determinants of their expression, involving molecules related to exocytosis and Rho signaling. Next, we wanted to examine the generality of our central findings, including that of (i) the identity of ubiquitously expressed CAMs and of (ii) the existence and identity of synaptic early- (subgraph 1) and late-gene interaction networks (subgraph 2). To this end, we performed single-cell RNAseq experiments to analyze electrophysiologically defined BS-PYR and RS-PYR neurons in the subiculum (Fig. 7A), which are the major projection targets of CA1-PYR cells. We hypothesized that their gene expression profiles are different, because they exhibit distinct excitability and synaptic properties (12, 13, 15).

We detected 9,671 ± 278 and 9,499 ± 405 genes in BS-PYR and RS-PYR cells, respectively (P = 0.85; Fig. 7B and Fig. S7A), with more consistent gene expression patterns than among CA1-PYRs (Fig. 7C and Fig. S7B). Next, we examined the molecular profiles of these cells, especially seeking exclusive markers that would unequivocally identify these functionally different cell types. However, we found that they were largely identical and detected only a surprisingly small number of exclusively expressed genes (Fig. 7D). Such genes included, for example, the neuropeptide cortistatin (Cort) and the c-Fos–induced growth factor (Figf) (a member of the VEGF family). Apart from these genes, BS-PYR and RS-PYR neurons expressed similar patterns of voltage-gated ion channels (Fig. 7E; see differently expressed ion channels in Fig. S7C), CAMs (Fig. 7F), ligand-gated ion channels, and exocytosis- and RhoGAP-/RhoGEF-related molecules (Fig. S7 C–G and Dataset S1).

We next examined expression of ubiquitously present CAMs (as identified in CA1 interneurons and PYR cells), which we reliably detected in subiculum PYR cells similar to CA1-PYRs (Figs. 4C and 7G), further strengthening the notion of general importance of these genes in neuron and synapse function. Although 42 CAMs had different expression levels between the two subiculum PYR cell types (P < 0.05; pairwise MW test; Dataset S1), none of these CAMs was exclusively expressed in one or the other cell type. However, their expression differed markedly from that of cell type-specifically expressed CAMs in CA1 PYR cells. Subiculum PYR cells consistently used CAMs that were not observed in CA1-PYR cells. For example, Clstn3 (significantly enriched only in interneurons in CA1), Ptpn5 and Ptprn2 (both were enriched in FS-PV cells), as well as Lrrc4 (not present in any CA1-PYR) were detected and highly expressed in subiculum PYR cells. Further cell-to-cell analysis corroborated molecular similarities between BS-PYR and RS-PYR cells, as confirmed by the inability of PCA (used on the complete transcriptome) to distinguish between these cells (Fig. 7H).

Finally, the subiculum neuron data allowed us to test the validity of synaptic gene correlation networks described in Fig. 6. We repeated the combined analysis of tissue and single-cell RNAseq data on subiculum PYR cells, independent from CA1 cells. In these analyses, we also identified two independent subgraphs, which were essentially identical to those found in CA1 RS, FS, and PYR cells. To reexamine the identity of the two networks, we singled out overlapping gene correlations in all five cell types (Fig. 7I) and found that they were nearly identical to subgraph 1 and subgraph 2, identified based on the three CA1 cell types (Fig. 6G). Therefore, these analyses independently confirmed the existence of coregulated and coexpressed synaptic gene networks, which may be generally present in cells in the brain independent of cell type identity.

All animal protocols and husbandry practices were approved by the Institutional Animal Care and Use Committee at Stanford University. Hippocampal slices (300 µm) were prepared from 3- to 4-wk-old wild-type CD1 mice, as described in ref. 10 and SI Materials and Methods.

Single-cell mRNA was performed using the Clontech’s SMARTer Ultra Low Input RNA Kit. As a first step, cells were collected via pipette aspiration into sample collection buffer, were spun briefly, and were snap frozen on dry ice. Samples were stored at −80 °C until further processing, which was performed according to manufacturer’s protocol. Library preparation was performed using Nextera XT DNA Sample Preparation Kit (Illumina) as described in the protocol. Then, cells were pooled and sequenced in an Illumina NextSeq500 instrument using 2 × 75 paired end reads on a NextSeq high-output kit (Illumina; see SI Materials and Methods for details).

After de-multiplexing the raw reads to single-cell datasets, we used Prinseq to remove short reads. After trimming and removal of overrepresented sequences and adapters, remaining reads were aligned to the mm10 genome with STAR aligner. Aligned read were converted to gene counts using HTSeq (see SI Materials and Methods for detailed parameter descriptions).

For gene expression analysis, we examined six functionally related categories: CAMs, voltage-gated ion channels, ligand-gated ion channels, synaptic exocytosis-related molecules, as well as RhoGAP and RhoGEF signaling-related molecules. For each of these, we assembled a comprehensive list, which included all genes examined (see SI Materials and Methods for a detailed description of categories).

All data analyses were performed using Mathematica10 (Wolfram Research). These analyses included (i) normalization of gene expression data, (ii) quality control, (iii) analysis of physiological properties, (iv) analysis of exon junctions, and (v) correlation and graph analyses of expression data (see SI Materials and Methods for more details).

Hippocampal slices (300 µm) were prepared from 3- to 4-wk-old wild-type CD1 mice incubated at 33 °C in sucrose-containing artificial cerebrospinal fluid (sucrose-ACSF) (85 mM NaCl, 75 mM sucrose, 2.5 mM KCl, 25 mM glucose, 1.25 mM NaH₂PO₄, 4 mM MgCl₂, 0.5 mM CaCl₂, and 24 mM NaHCO₃) for 1 h, and then held at room temperature until recording. Cells were visualized by infrared differential interference contrast optics in an upright microscope (Olympus; BX-61WI). Recordings were performed using borosilicate glass pipettes with filament (Sutter Instruments; BF100-58-10; o.d., 1.0 mm; i.d., 0.58 mm; 10-cm length) at 33 °C in ACSF (126 mM NaCl, 2.5 mM KCl, 10 mM glucose, 1.25 mM NaH₂PO₄, 2 mM MgCl₂, 2 mM CaCl₂, and 26 mM NaHCO₃) with a standard intracellular solution (95 mM K-gluconate, 50 mM KCl, 10 mM Hepes, 4 mM Mg-ATP, 0.3 Na-GTP, 10 mM phosphocreatine; pH 7.2, 300 mOsm).

Single-cell mRNA was performed using the Clontech’s SMARTer Ultra Low Input RNA Kit for Sequencing, version 3 (for tissue/RS-INT/FS-INT/CA1-PYR samples) and version 4 (for BS-PYR/RS-PYR samples). As first step, cells were collected via pipette aspiration into 1.1 µL of 10× collection buffer, were spun briefly, and were snap frozen on dry ice. Samples were stored at −80 °C until further processing, which was performed according to manufacturer’s protocol. Resulting cDNA was harvested and analyzed on the Fragment analyzer automated fragment analyzer by Advanced Analytical. Only cells that had a concentration higher than 0.05 ng/µL were selected for library preparation. Library preparation was performed using Nextera XT DNA Sample Preparation Kit (Illumina) as described in the protocol. Following library preparation, cells were pooled and sequenced in an Illumina NextSeq500 instrument using 2 × 75 paired end reads on a NextSeq high-output kit (Illumina).

After demultiplexing the raw reads to single-cell datasets, we used Prinseq to remove short reads (-min_len 30), to trim the first 10 bp on the 5′-end (-trim_left 10), to trim reads with low quality on the 3′-end (-trim_qual_right 25), and to filter low complexity reads (-lc_method entropy \-lc_threshold 65). We used FASTQC to determine overrepresented sequences and removed those using cutadapt (-e 0.15 –m 30). We then used Prinseq to remove orphan pairs less than 30 bp in length followed by removal of Nextera adapters using Trim Galore (–stringency 1). Remaining reads were aligned to the mm10 genome with STAR using the following options (–outFilterType BySJout \–outFilterMultimapNmax 20 \–alignSJoverhangMin 8 \–alignSJDBoverhangMin 1 \–outFilterMismatchNmax 999 \–outFilterMismatchNoverLmax 0.04 \–alignIntronMin 20 \–alignIntronMax 1000000 \–alignMatesGapMax 1000000 \–outSAMstrandField intronMotif). Then, in each sample, aligned reads were converted to counts for every gene using HTSeq (-m intersection-nonempty \-s no).

For gene expression analysis, we examined six functionally related categories: CAMs, voltage-gated ion channels, ligand-gated ion channels, synaptic exocytosis-related molecules, as well as RhoGAP and RhoGEF signaling-related molecules. For each of these, we assembled a comprehensive list, which included all genes examined (see Dataset S1 for more details).

All data analyses were performed using Mathematica10 (Wolfram Research). These analyses included (i) normalization of gene expression data, (ii) quality control (QC), (iii) analysis of physiological properties, (iv) analysis of exon junctions, and (v) correlation and graph analyses of expression data.