Research Article

Three-dimensional intact-tissue sequencing of single-cell transcriptional states

See allHide authors and affiliations

Science  27 Jul 2018:
Vol. 361, Issue 6400, eaat5691
DOI: 10.1126/science.aat5691

Transcriptome mapping in the 3D brain

RNA sequencing samples the entire transcriptome but lacks anatomical information. In situ hybridization, on the other hand, can only profile a small number of transcripts. In situ sequencing technologies address these shortcomings but face a challenge in dense, complex tissue environments. Wang et al. combined an efficient sequencing approach with hydrogel-tissue chemistry to develop a multidisciplinary technology for three-dimensional (3D) intact-tissue RNA sequencing (see the Perspective by Knöpfel). More than 1000 genes were simultaneously mapped in sections of mouse brain at single-cell resolution to define cell types and circuit states and to reveal cell organization principles.

Science, this issue p. eaat5691; see also p. 328

Structured Abstract

INTRODUCTION

Single-cell RNA sequencing has demonstrated that both stable cell types and transient cell states can be discovered and defined by transcriptomes. In situ transcriptomic methods can map both RNA quantity and position; however, it remains challenging to simultaneously satisfy key technological requirements such as efficiency, signal intensity, accuracy, scalability to large gene numbers, and applicability to three-dimensional (3D) volumes. Well-established single-molecule fluorescence in situ hybridization (FISH) approaches (such as MERFISH and seqFISH) have high detection efficiency but require long RNA species (more than 1000 nucelotides) and yield lower intensity than that of enzymatic amplification methods (tens versus thousands of fluorophores per RNA molecule). Other pioneering in situ sequencing methods (via padlock probes and fluorescent in situ sequencing) use enzymatic amplification, thus achieving high intensity but with room to improve on efficiency.

RATIONALE

We have developed, validated, and applied STARmap (spatially-resolved transcript amplicon readout mapping). STARmap begins with labeling of cellular RNAs by pairs of DNA probes followed by enzymatic amplification so as to produce a DNA nanoball (amplicon), which eliminates background caused by mislabeling of single probes. Tissue can then be transformed into a 3D hydrogel DNA chip by anchoring DNA amplicons via an in situ–synthesized polymer network and removing proteins and lipids. This form of hydrogel-tissue chemistry replots amplicons onto an optically transparent hydrogel coordinate system; then, to identify and quantify RNA species-abundance manifested by DNA amplicons, the identity of each species is encoded as a five-base barcode and read out by means of an in situ sequencing method that decodes DNA sequence in multicolor fluorescence. Using a new two-base sequencing scheme (SEDAL), STARmap was found to simultaneously detect more than 1000 genes over six imaging cycles, in which sequencing errors in any cycle cause misdecoding and are effectively rejected.

RESULTS

We began by (i) detecting and quantifying a focused 160-gene set (including cell type markers and activity-regulated genes) simultaneously in mouse primary visual cortex; (ii) clustering resulting per-cell gene expression patterns into a dozen distinct inhibitory, excitatory, and non-neuronal cell types; and (iii) mapping the spatial distribution of all of these cell types across layers of cortex. For validation, per-cell-type gene expression was found to correlate well both with in situ hybridization results and with single-cell RNA sequencing, and widespread up-regulation of activity-regulated genes was observed in response to visual stimulation. We next applied STARmap to a higher cognitive area (the medial prefrontal cortex) and discovered a more complex distribution of cell types. Last, we extended STARmap to much larger numbers of genes and spatial scales; we measured 1020 genes simultaneously in sections—obtaining results concordant with the 160-gene set—and measured 28 genes across millimeter-scale volumes encompassing ~30,000 cells, revealing 3D patterning principles that jointly characterize a broad and diverse spectrum of cell types.

CONCLUSION

STARmap combines hydrogel-tissue chemistry and in situ DNA sequencing to achieve intact-tissue single-cell measurement of expression of more than a thousand genes. In the future, combining this intact-system gene expression measurement with complementary cellular-resolution methodologies (with which STARmap is designed to be compatible)—including in vivo activity recording, optogenetic causal tests, and anatomical connectivity in the same cells—will help bridge molecular, cellular, and circuit scales of neuroscience.

STARmap for 3D transcriptome imaging and molecular cell typing.

STARmap is an in situ RNA-sequencing technology that transforms intact tissue into a 3D hydrogel-tissue hybrid and measures spatially resolved single-cell transcriptomes in situ. Error- and background-reduction mechanisms are implemented at multiple layers, enabling precise RNA quantification, spatially resolved cell typing, scalability to large gene numbers, and 3D mapping of tissue architecture.

Abstract

Retrieving high-content gene-expression information while retaining three-dimensional (3D) positional anatomy at cellular resolution has been difficult, limiting integrative understanding of structure and function in complex biological tissues. We developed and applied a technology for 3D intact-tissue RNA sequencing, termed STARmap (spatially-resolved transcript amplicon readout mapping), which integrates hydrogel-tissue chemistry, targeted signal amplification, and in situ sequencing. The capabilities of STARmap were tested by mapping 160 to 1020 genes simultaneously in sections of mouse brain at single-cell resolution with high efficiency, accuracy, and reproducibility. Moving to thick tissue blocks, we observed a molecularly defined gradient distribution of excitatory-neuron subtypes across cubic millimeter–scale volumes (>30,000 cells) and a short-range 3D self-clustering in many inhibitory-neuron subtypes that could be identified and described with 3D STARmap.

In biological tissues, diversity of function arises from diversity of form—in part, via the complexity of cell-specific gene expression, which defines the distinct three-dimensional (3D) molecular anatomy and cellular properties of each tissue. In situ transcriptomic tools for the spatial mapping of gene expression with subcellular resolution have emerged that may be applicable to probing these tissue structure–function relationships, including both multiplexed in situ RNA hybridization and in situ RNA sequencing (110). Current in situ sequencing approaches face the challenge of implementing enzymatic reactions in the dense, complex tissue environment and currently suffer from low efficiency (2), but the potential value of such intact-tissue sequencing could be enormous; in comparison with hybridization-based multiplexing/readout, which uses multiple polynucleotide probes to encode gene identity (35), sequencing operates with single-nucleotide resolution and thus inherently provides greater information. In addition, in situ sequencing methods typically use signal amplification, which is important for the detection of short transcripts (such as neuropeptides) and for high-quality imaging in thick tissue blocks. However, current sequencing methods have not yet been successfully applied to 3D volumes of intact tissue because of fundamental limitations in requisite sensitivity, fidelity, and scalability for throughput in tissues such as the mammalian brain.

Hydrogels have been widely used for extracellular 3D scaffolding in applications across biology and medicine (1113). Recently developed hydrogel-tissue chemistry (HTC) methodologies (14), beginning with CLARITY (15), physically link in situ–synthesized polymers with selected intracellular biomolecules. This process transforms the tissue, from within its constituent cells, into a new state suitable for high-resolution volumetric imaging and analysis compatible with many kinds of molecular phenotyping for proteins, nucleic acids, and other targets (15). HTC-based hydrogel-embedding strategies have been extended to nucleic acid analyses in the form of in situ hybridization for RNA (1619), but these have not yet been extended to in situ RNA sequencing—which would have the potential to reveal the full molecular complexity of the transcriptome. In nontissue environments, however, purely synthetic hydrogels have been used to accommodate enzymatic reactions that include DNA sequencing (20), and if biological tissue could be converted into a hydrogel-embedded form compatible with creation, retention, and functional presentation of RNA-derived or hybridized complementary DNA (cDNA), it might be possible to perform 3D in situ sequencing within such a tissue-hydrogel formulation—leveraging the crucial attendant properties of optical transparency, reduced background, elevated diffusion rate, and greater mechanical stability. We achieved this goal with the development and application of a sequencing-based method (spatially-resolved transcript amplicon readout mapping, or STARmap) for targeted 3D in situ transcriptomics in intact tissue (Fig. 1A); using STARmap, we were able to identify organizational principles of a full spectrum of cell types, which would not have been otherwise accessible for identification in the adult mammalian brain.

Fig. 1 STARmap principles: in situ RNA sequencing for spatial transcriptomics within the 3D tissue environment.

(A) STARmap overview schematic. After brain tissue is prepared (mouse brain protocols are available in the supplementary materials, materials and methods), the custom SNAIL probes that encounter and hybridize to intracellular mRNAs (dashed lines) within the intact tissue are enzymatically replicated as cDNA amplicons. The amplicons are constructed in situ with an acrylic acid N-hydroxysuccinimide moiety modification (blue) and then copolymerized with acrylamide to embed within a hydrogel network (blue wavy lines), followed by clearance of unbound lipids and proteins (fig. S2). Each SNAIL probe contains a gene-specific identifier segment (red) that is read-out through in situ sequencing with two-base encoding for error correction (SEDAL) (fig. S3). Last, highly multiplexed RNA quantification in three dimensions reveals gene expression and cell types in space. (B) SNAIL logic. A pair of primer and padlock probes amplifies target-specific signals and excludes noise known to commonly arise from nonspecific hybridization of a single probe. (C and D) Only adjacent binding of primer and padlock probes leads to signal amplification. mRNA A represents Gapdh, and mRNA B represents Actb. Both fluorescent images show Gapdh (gray) mRNA and cell nuclei (blue) labeling in mouse brain slice; there is an absence of labeling with mismatched primer and padlock (right). Scale bar, 10 μm. (E) In situ sequencing of DNA amplicons in the tissue-hydrogel complex via SEDAL, the sequencing-by-ligation method devised for STARmap. For each cycle, the reading probes (gray line without star-symbol label) contain an incrementally increasing-length run of degenerate bases (N representing an equal mixture of A, T, C, and G) with phosphate at the 5′ end (5′P) to set the reading position; the decoding probes (gray line with star-symbol label) are labeled by fluorophores with color coding for the dinucleotide at the 3′ end. Only if both probes are perfectly complementary to the DNA template (black lower sequence) can the two kinds of probes then be ligated to form a stable product with a high melting temperature, allowing later imaging after unligated probes are washed away. After each imaging cycle, probes are stripped away from the robust tissue-hydrogel by using 60% formamide so that the next cycle can begin. X, unknown base to be read; red underline, decoded sequence; Ch1 to Ch4, fluorescence channels. Scale bar, 2 μm.

Results

Design and validation of STARmap principles

One component is an efficient approach for in situ amplification of a library of cDNA probes hybridized with cellular RNAs (this approach is termed SNAIL, for specific amplification of nucleic acids via intramolecular ligation). Reverse transcription may be the major efficiency-limiting step for in situ sequencing (7, 21), and SNAIL bypasses this step with a pair of primer and padlock probes (fig. S1A) designed so that only when both probes hybridize to the same RNA molecule can the padlock probe be circularized and rolling-circle-amplified to generate a DNA nanoball (amplicon) that contains multiple copies of the cDNA probes (Fig. 1, A to D). This mechanism ensures target-specific signal amplification and excludes noise that invariably otherwise arises from nonspecific hybridization of single probes. The outcome includes much higher absolute intensity and signal-to-noise ratio (SNR) as compared with those of commercial single-molecule fluorescent in situ hybridization (smFISH) probes (fig. S1, B to F) and substantial improvement of detection efficiency (comparable with that of single-cell RNA sequencing), with simplified experimental procedures compared with previous in situ RNA sequencing methods (fig. S1, G to I).

To enable cDNA amplicon embedding in the tissue-hydrogel setting, amine-modified nucleotides were spiked into the rolling-circle amplification reaction, functionalized with an acrylamide moiety by using acrylic acid N-hydroxysuccinimide esters, and copolymerized with acrylamide monomers so as to form a distinct kind of hydrogel-DNA amplicon network (Fig. 1A and fig. S2A). The resulting tissue-hydrogel was then subjected to protein digestion and lipid removal in order to enhance transparency (fig. S2, B to E). This design chemistry dictates that amplicons are covalently linked with the hydrogel network, and such cross-linking is essential to maintain the position and integrity of the amplicons through many cycles of detection (fig. S2, F to H).

A five-base barcode (library size of 1024) was designed and built into each padlock probe as a gene-unique identifier to be sequenced, thus enabling multiplexed gene detection (Fig. 1A). Sequencing-by-synthesis paradigms were avoided because these require elevated reaction temperatures, which in turn are problematic for high-resolution imaging and sample stability (16) in comparison with sequencing-by-ligation methods that can be implemented at room temperature. However, none of the reported or commercially available sequencing-by-ligation methods exhibit the necessary SNR or accuracy for this challenging intact-tissue application: Supported Oligo Ligation Detection (SOLiD) sequencing causes strong background fluorescence in biological samples (10), whereas combinatorial probe anchor ligation (cPAL) sequencing (22) lacks an error-rejection mechanism (fig. S3). For this reason, an approach we term sequencing with error-reduction by dynamic annealing and ligation (SEDAL) was devised specifically for STARmap (fig. S3).

SEDAL uses two kinds of short, degenerate probes: reading probes to decode bases, and fluorescence probes to transduce decoded sequence information into fluorescence signals. The two short probes only transiently bind to the target DNA and ligate to form a stable product for imaging only when a perfect match occurs; after each cycle corresponding to a base readout, the fluorescent products are stripped with formamide, which eliminates error accumulation as sequencing proceeds (Fig. 1E and fig. S3B). In contrast to SOLiD, SEDAL exhibits minimal background (fig. S3, C to F). A two-base encoding scheme was designed and implemented in order to mitigate any residual errors related to imaging high densities of spots (fig. S3, G and H). On the basis of a panel of four very highly expressed test genes in mouse brain (to mimic amplicon crowdedness as would be encountered in highly multiplexed gene-detection), we found that the error rate of STARmap was more than an order of magnitude lower than prior cPAL methods (~1.8 versus 29.4%) (fig. S3, I to L) (17).

Spatial cell typing in primary visual cortex with 160-gene STARmapping

To test whether STARmap could deliver on the initial goal of high-content 3D intact-tissue sequencing of single-cell transcriptional states with the necessary sensitivity and accuracy, we applied STARmap to a pressing challenge in neuroscience: detecting and classifying cell types and corresponding tissue-organization principles in the neocortex of the adult mouse brain. The anatomy and function of the mouse primary visual neocortex have been extensively studied (23), a setting which here allows validation of our results by comparison with prior findings that span multiple papers, methodologies, and data sources (but the full diversity of deeply molecularly defined cell types within the visual cortex has not yet been spatially resolved in a single experiment, precluding identification of potentially fundamental joint statistics and organizational principles across 3D volumes). Among many examples of the experimental leverage such information could provide, joint 3D cell-typology mapping might be used to help decode the spatiotemporal logic of neural activity–triggered gene expression as a function of cell type and spatial location.

We therefore used five-base barcoded SNAIL probes over six rounds of in situ SEDAL sequencing in coronal mouse brain slices (Figs. 1A and 2, A and B) to survey a large but focused and curated gene set [160 genes including 112 putative cell-type markers collated from mouse cortical single-cell RNA sequencing (24, 25) and 48 activity-regulated genes (ARGs) (26, 27)]. In one arm of the experiment, visually evoked neural activity was provided to a cohort of mice via 1 hour of light exposure after 4 days of housing in the dark; other mice were kept continuously in the dark (27, 28). Eight-μm-thick volumes containing up to 1000 cells covering all cortical layers were imaged. After six rounds of sequencing, fluorescent Nissl staining was used to segment cell bodies, allowing attribution of amplicons to individual cells (fig. S4, A and B). The values corresponding to amplicons-per-cell and genes-per-cell varied substantially (Fig. 2C), whereas the 160-gene expression pattern was consistent between biological replicates [correlation coefficient (r) = 0.94 to 0.95] (Fig. 2D), revealing reliable detection of transcript diversity at the single-cell level. Because only 160 genes were encoded out of the 1024 possible barcodes from five bases, we were able to quantify sequencing errors that resulted in sequences being corrupted from the 160 true barcodes to the 864 invalid barcodes, which was remarkably low at 1 to 4%. We found that this 160-gene pilot faithfully reproduced the spatial distribution of known cortical layer markers and interneurons, illustrated here via comparison of in situ images from paired public atlases (29) and STARmap results (Fig. 2E).

Fig. 2 STARmapping cell types in V1.

(A) Experimental design. Mice were dark housed, before sacrifice, for 4 days and then either kept in the dark or exposed to light for 1 hour. V1 was coronally sectioned, and RNAs of 112 cell type markers and 48 activity-regulated genes were quantified by means of STARmap. (B) Raw fluorescence images of in-process STARmap with the full view of cycle 1 (top) and zoomed views across all six cycles (bottom). Full field: 1.4 by 0.3 mm, scale bar, 100 μm; zoomed region: 11.78 by 11.78 μm, scale bar, 2 μm; Channel, color code for the four fluorescence channels; L1 to L6, the six neocortical layers; cc, corpus callosum; HPC, hippocampus. (C) Histograms. Shown are detected reads (DNA amplicons) per cell (left), and genes per cell (right). (D) Quantitative reproducibility of biological replicates, whether in the light or dark condition: log2(amplicon quantity) for 160 genes across the whole imaging region plotted. Rep1, expression value in first replicate; rep2, expression value in second replicate. (E) Validation of STARmap. (Left) in situ images from Allen Institute of Brain Science (AIBS). (Right) RNA pattern of individual genes extracted from 160-gene STARmap, which reliably reproduced the spatial gene expression pattern from AIBS. (F) Uniform manifold approximation plot (UMAP), a nonlinear dimensionality reduction technique used to visualize the similarity of cell transcriptomes in two dimensions, showing consistent clustering of major cell types across 3142 cells pooled from four biological replicates: 2199 excitatory neurons, 324 inhibitory neurons, and 619 non-neuronal cells. (G) Gene expression heatmap for 112 cell-type markers aligned with each cell cluster, showing clustering by inhibitory, excitatory, or non-neuronal cell types. Expression for each gene is z-scored across all genes in each cell. (H) Representative cell-resolved spatial map in neocortex and beyond. Cell types are color-coded as in (F). (I to N) Clustering of excitatory and inhibitory subtypes. [(I) and (L)] UMAP plots, [(J) and (M)] bar plots of representative genes (mean ± 95% confidence interval expression across all cells in that cluster, with each bar scaled to the maximum mean expression across all clusters), and [(K) and (N)] in situ spatial distribution of [(I) to (K)] excitatory and [(L) to (N)] inhibitory neurons. The number of cells in each cluster was as follows: eL2/3, 589; eL4, 649; eL5, 393; eL6, 368; PV neurons, 111; VIP neurons, 46; SST neurons, 46; and NPY neurons, 56. Inclusion of cells in clusters was guided entirely by amplicon representation in each cell without using spatial information; excitatory cell clusters were then named according to the spatial layering observed for that cluster, whereas inhibitory cell clusters were named according to the dominant cell-type amplicon based on the strong segregation of amplicon markers.

We next performed cell classification using expression data of the 112 cell-type markers. First, >3000 cells pooled from four biological replicates were clustered into three major cell types (excitatory neurons, inhibitory neurons, and non-neuronal cells) by using graph-based clustering after principal-component decomposition (30) and then further subclustered under each category (Fig. 2, F to H, and fig. S4C). The richly defined excitatory neurons segregated into four major types (here denoted eL2/3, eL4, eL5, and eL6) by spatial correspondence with anatomic cortical layers and expression profiles of known layer-specific gene markers (Fig. 2, I to K, and fig. S5, A and B). Although spatial organization of the four excitatory types exhibited a layered pattern, there was extensive intermixing among different cell types within each layer. Inhibitory neurons were also clustered into four major types, here denoted by the dominant interneuron marker of each subtype [VIP, SST, NPY, and PV (Pvalb)] (Fig. 2, L to N, and fig. S5, C and D); the VIP and NPY type were observed to distribute more to the upper layers (L1 to L3), whereas SST and PV types were found more commonly in the lower layers (L4 to L6). We also detected non-neuronal cell types, including astrocytes, oligodendrocytes, endothelial cells, and smooth muscle cells (fig. S6). The number of major cell types illustrated here (12 in total) can be further broken down (single-cell RNA sequencing can lead to classification into 40 or more subtypes, which is consistent with the readily apparent heterogeneity of gene expression within each type) (figs. S5 and S6). With our targeted 112-gene set and at the size of 600 to 800 cells per sample, all 12 major cell types could be reliably detected without batch effects with highly similar spatial patterning among four biological replicates (defined as samples prepared from different animals) (fig. S7) and matched with published single-cell RNA sequencing results (fig. S8).

We next sought to take advantage of STARmap’s quantitative capabilities at the single-cell level, in order to test differential gene expression analyses across experimental conditions, in molecularly defined cell types. To this end, we assessed visual stimulus–dependent gene expression patterns (via 48 defined ARGs with single-cell resolution in situ). Further developing the single-cell RNA sequencing procedure, mouse brains were flash-frozen with minimal handling time after sacrifice (<5 min), for maximal preservation of native transcriptional signatures. Consistent with prior reports (2628), we observed global induction of known immediate-early genes (Fos, Egr1, and Egr2) (Fig. 3A) in the primary visual cortex (V1) upon 1 hour of light exposure. At single-cell resolution, the quantitative extent (fold change in expression) of ARG changes exhibited striking diversity across neuronal cell types (Fig. 3, B and C, and fig. S9) (28). In general, ARG expression programs in excitatory neurons across different layers were highly similar, whereas ARG expression programs in inhibitory cells exhibited much more distinct cell type–specific characteristics (fig. S9C); for example, Egr2 exhibited light-induction across excitatory neurons (Fig. 3D) but not in inhibitory neurons, whereas Prok2 was up-regulated in Vip inhibitory neurons (Fig. 3C) (22). Last, because neural activity can trigger cotranscription of noncoding RNAs from within enhancers of ARGs (26, 31), we also studied exemplars of these enhancer RNAs (eRNAs) (here, eRNAs 1 to 5 of the Fos gene); these transcripts, not polyadenylated, would be very difficult to measure with current single-cell RNA sequencing. eRNA3 was identified as the most notable and consistent ARG marker (fig. S9B).

Fig. 3 STARmapping behavioral experience: Detecting and quantifying cell type–specific regulation of ARGs.

(A) Validation. Shown is spatial expression pattern in the visual cortex of prototypical ARGs known as immediate early genes. Sacrifice was in darkness or after 1 hour of light exposure. (B and C) Volcano plots of log fold-change in gene expression between light and dark conditions in inhibitory and excitatory cell types. Genes with significantly increased or decreased expression (false discovery rate–adjusted P < 0.05, Wilcoxon rank-sum test) are labeled in green, and the most significantly changed genes (P < 0.05 and fold change > 2) are labeled in red. Many ARGs showed cell-type specificity, pointing to discovery of unanticipated cell type–specific logic of excitation-transcription coupling. (D) Violin plot of Egr2 expression by cell type. ****P < 0.0001, n.s. not significant, Wilcoxon rank-sum test; red-labeled cell types, fold change >2.

Comparing spatial cell-type distributions in frontal and sensory cortices

We then investigated to what extent the cell types of the higher cognitive cortex resemble those of the sensory cortex, as exemplified by V1. We applied the same 160-gene set to STARmapping the cell populations of the medial prefrontal cortex (mPFC) (Fig. 4A), which is involved in high-level cognitive functions such as attention and memory and is thought to be dysregulated in major psychiatric disorders (32). We identified 15 distinct molecular cell types, including six excitatory neuron subtypes (eL2/3, eL5-1, eL5-2, eL5-3, eL6-1, and eL6-2, annotated by anatomic cortical layers), five inhibitory neuron subtypes (VIP, Reln, SST, Lhx6 and NPY, annotated by dominant gene markers), and four non-neuronal types (astrocytes, oligodendrocytes, endothelial cells, and smooth-muscle cells) (Fig. 4B and fig. S10).

Fig. 4 STARmapping cell types and neural activity in mPFC.

(A) Diagram of targeted region (red box) containing primarily prelimbic cortex (PrL) within mPFC. (B) UMAP visualization of all inhibitory (VIP, Reln, SST, Lhx6, and NPY), excitatory (eL2/3, eL5-1, eL5-2, eL5-3, eL6-1, and eL6-2), and non-neuronal (Astro, Oligo, Smc, and Endo) cell types. (C) Spatial visualization of cell type layout in mPFC, using the same color scheme as in (B). (D) Barplot of representative genes per cluster (mean ± 95% confidence interval), with each bar scaled to the maximum mean expression for that gene across clusters. (E) Piecharts showing the relative proportion of each major and minor cell type in both mPFC and visual cortex. (F) Violin plots of Fos gene induction in different excitatory cell types in mPFC in response to cocaine. The mice were sacrificed after 1 hour of cocaine or saline injection. Expr, normalized expression; n.s., not significant; *P < 0.05, ***P < 0.001, ****P < 0.0001, likelihood ratio test. Astro, astrocytes; Oligo, oligodendrocytes; Smc, smooth muscle cells; Endo, endothelial cells.

The spatial organization of broad cell types in mPFC resembled that of V1 with intermixed excitatory neuronal layers and sparsely distributed inhibitory neurons (Fig. 4C); however, the nature and composition of neuronal subtypes in mPFC and V1 strikingly differed (Fig. 4, D and E). For excitatory subtypes, mPFC lacks eL4 (which is consistent with previous reports) (33) and exhibits reduced eL2/3 and vast expansion of eL5 and eL6 compared with that of V1 (Fig. 4E). Many new types of cell were discovered, including three eL5 subtypes and two eL6 subtypes, as characterized by gene markers Sema3e, Plcxd2, Tpbg, Syt6, and Ctgf, respectively (Fig. 4D).

Substantially different tissue organization by cell type was also observed for inhibitory subclusters. Sst-, Vip-, and Npy-positive subtypes in mPFC were represented similarly among all inhibitory neurons compared with those in V1, whereas Pvalb-positive cells were comparatively much sparser. In V1, Reln-positive neurons coexist with Sst and Npy, whereas in mPFC, these segregate as a single cluster, with ~50% comarked by Ndnf; we also discovered a new inhibitory subtype labeled by Lhx6, which in fact constitutes the most abundant inhibitory subtype in mPFC (Fig. 4E). Although the 5-HT(3A) receptor (Htr3a) expression has been reported in cortical inhibitory neurons (34), Htr3a has not been ranked as a critical genetic marker of inhibitory subtypes in V1. In mPFC, however, we found that Htr3a distinguishably marks a large fraction of Vip+ neurons and a subset of Reln+ neurons (fig. S10D).

Superficial layers (L1 to L3) were found to contain Vip, Reln, and Npy subtypes, whereas deeper layers (L5 to L6) were found to contain all of the inhibitory subtypes. All of the 15 cell types with tissue-level spatial organization could be reliably detected with STARmap across four biological replicates (fig. S11). The capability of STARmap for multidimensional cell typing in mPFC was further demonstrated in the setting of activity dependence, supporting the possibility of defining cell types in part by communication properties, including activity during behavior (35, 36). One hour after cocaine injection (37), a specific subpopulation of deep-layer excitatory neurons (such as Tpbg labeled eL5-2) in mPFC was activated compared with that in saline-injected control mice (Fig. 4F), revealing STARmap capability for identifying functional segregation of neuronal subtypes in mPFC.

Scaling STARmap to more than 1000 genes

To further test the scalability of STARmap, we extended our gene list from 160 to 1020 genes, leveraging previously published single-cell RNA sequencing data (24). The 1020-gene set was first validated in mouse hippocampal neuron culture, with successful resolution of neuronal and glial cells (fig. S12). We then probed mouse V1 neocortex with the 1020-gene set in order to evaluate performance in spatial cell typing in comparison with that of the 160-gene set. Amplicons obtained in the 1020-gene experiment were much denser in cells as compared with those in 160-gene experiments but were optically resolvable in 3D with high-resolution imaging and postimaging deconvolution (Fig. 5A).

Fig. 5 Simultaneous mapping of 1020 genes in V1 by STARmap.

(A) Input fluorescence data. (Left) Maximum-intensity projection of the first sequencing round for 1020 gene experiment, showing all four channels simultaneously. Yellow square, zoom region. Scale bar, 100 μm. (Right) Zoom into a single cell showing spatial arrangement of amplicons in three dimensions across six sequencing rounds. (B) Joint UMAP plot showing all excitatory (HPC, eL2/3, eL4, eL5, eL6-1, and eL6-2), non-neuronal (Smc, Other, Olig, Micro, Endo, and Astro), and inhibitory (PVALB, SST, VIP, and NPY) cell types. (C) Plot of all differentially expressed genes across every cluster, with P < 10−12 and log fold change > 1.5. (D) Spatial map of all excitatory, non-neuronal, and inhibitory cell types in visual cortex using the same color code of (B). HPC, hippocampus; Smc, smooth muscle cells; Other, other unclassified cells; Oligo, oligodendrocytes; Micro, microglia; Endo, endothelia cells; Astro, astrocytes.

We observed that a higher percentage (40%) of amplicons were filtered out in the 1020-gene experiments by our error-rejection mechanism (fig. S3H) in comparison with the four-gene experiments (20%) (fig. S3L), indicating that a more frequent initial color-misassignment potentially resulted from amplicon merging or optical resolution and further demonstrating the importance of our designed error-rejection mechanism. Crucially, despite the read loss, we successfully clustered single cells of the imaging area into 15 annotated cell types and one unclassified type using 1020 genes and the same data analysis pipeline from the focused 160 gene probe set (Fig. 5, B and C, and fig. S13). Three new cell types were identified in addition to the 12 cell types detected by 160 genes (Fig. 5B): eL6 was resolved into two subtypes (eL6-1 and eL6-2), a previously uncharacterized hippocampal excitatory subtype (HPC) was identified, and microglial cells were cleanly identified with an expansion of non-neuronal cell type markers in the 1020-gene set.

Beyond those advances, the 1020-gene findings also successfully reproduced the cell types (and their spatial patterning) from the 160-gene findings and further allowed discovery of multiple new gene markers for each cell type (for example, 3110035E14Rik for deep layers, Cnot6l for Sst neurons, and Cplx1 for Pvalb neurons) (Fig. 5D and fig. S13). These molecularly defined cell types were highly reproducible between biological replicates for 1020-gene detection and were concordant with published single-cell RNA sequencing results (fig. S14). We further assessed the possibility of scaling up STARmap to accommodate higher gene numbers; although the STARmap scheme can encode and decode more than 1 million codes and the physical volume of mammalian cells is not limiting for amplification of more than 1000 genes (fig. S15), the 1020-gene experiments approached the upper limit of the optical volume of cells (fig. S15E); for those cases in which more genes are needed, STARmap may cover the whole transcriptome with more sequencing rounds of serial 1000-gene detection, or via optical resolution enhanced with super-resolution microscopy (38, 39) or the physical swelling typical of the hydrogel-tissue chemistries (14, 19).

Adapting STARmap to thick tissue blocks for 3D analyses

In neuroscience, addressing the 3D complexity of both neurons and neural circuits has generally required the development and use of thick tissue blocks or fully intact brains for functional and structural readouts, including electrophysiology, imaging of activity, and analysis of morphology and connectivity. Therefore, for linking these readout measures from intact or semi-intact tissue preparations with cellular-resolution gene expression readouts from the very same preparations, methods of 3D spatial transcriptomic analysis in thick tissues have long been sought in order to achieve datastream registration as well as preserve 3D morphology and to obtain readouts from very much larger cell numbers (2). The initial experiments were carried out in brain slices no more than one cell body thick; we therefore next developed and tested STARmap to overcome limitations in diffusional access and imaging throughput for large tissue volumes, with a modified strategy for linearly reading out gene expression at cellular resolution so as to enable high-throughput molecular analysis in tissue volumes (Fig. 6A and fig. S16). Specificity and penetration depth of large-volume STARmap were tested initially by using Thy1::YFP mouse brains, in which STARmap successfully detected yellow fluorescent protein (YFP) mRNA across 150 μm of tissue thickness and specifically colocalized YFP protein and mRNA at single-cell resolution (Fig. 6B) without labeling the tens of thousands of interspersed neighboring cells.

Fig. 6 3D architecture of cell types in visual cortex volumes.

(A) Volumetric STARmapping via sequential SEDAL gene readout. Using a modified STARmap procedure (fig. S16) and cyclic gene readout (four genes in each cycle), large tissue volumes can be rapidly mapped at single-cell resolution without oversampling each amplicon. (B) Validation showing specific STARMAP labeling of YFP-expressing neurons (from transgenic Thy1::YFP mouse line) in 3D cortical volume. Scale bar, 0.5 mm. (C) Representative labeling of (left) major cell types, (left center) layer-specific markers, (right center) inhibitory markers, and (right) activity-regulated genes acquired over multiple rounds in visual cortex STARmap volumes. (D) Per-cell expression matrix of 28 genes from 32,845 single cells from one volume clustered into multiple excitatory, inhibitory, and non-neuronal cell types, z-scored across genes for each cell in order to normalize for mean differences in total signal between cells. Columns are sorted by order of sequencing rounds as conducted, in groups of four. (E) (Top) Spatial histograms of excitatory, inhibitory, and non-neuronal cell types, using same color labels as (D). Cells were counted in 5-μm bins in a 2D max-projection and plotted in cell-count-per-micrometer units as a function of distance from the corpus callosum (cc) to the pia, averaged across the bins perpendicular to the cortical layers. (Bottom) Plot of max-projected cell locations color-coded by cluster as in (D). (F) Spatial distribution of each cell type (excitatory, inhibitory, and non-neuronal) and subtypes in three dimensions. Each dot represents a single cell; spatial dimensions are in micrometers. (G) Average nearest-neighbor distances computed in three dimensions between all excitatory cells (Excite) and each inhibitory cell type. For self-comparisons, the nearest neighbor was defined as the closest nonidentical cell; persistent self-correlation reveals self-clustering of inhibitory subtypes. (H) Same distances as (H) but using shuffled (randomized) cell-type labels. (I) Nearest-neighbor distances computed in three dimensions between each inhibitory cell of a certain type and any member of the same type (Inhib → Inhib, eg VIP → VIP) or any excitatory neuron (Inhib → Excite). **** P < 0.0001, Wilcoxon rank-sum test.

We then extended the spatial cell-typing of mouse V1 to more than 30,000 cells across volumes spanning all six layers and the corpus callosum. Using a curated gene set including 23 cell-type markers and five ARGs read out over seven cycles of linear SEDAL sequencing (Fig. 6, C and D, and fig. S17), we applied K-means clustering of marker genes (supplementary materials, materials and methods) for each cell type (recovering 11 cell types corresponding to the majority of those extracted by the 160-gene experiment). We found that 3D patterning of the 11 cell types (Fig. 6, E and F) was consistent with the 160-gene thin-section tissue findings but provided an accurate and quantitative profiling of cellular distribution across space, with much larger cell numbers. As reflected by both spatial-histogram (Fig. 6E) and correlational analyses (fig. S17B), excitatory subtypes exhibited a layered gradient distribution, with the spatial density of each subtype decaying across space into adjacent layers. By contrast, inhibitory subtypes were dispersed, albeit with layer preferences exhibited by the Vip subtype (largely located in layer 2/3) and the Sst and Pvalb subtypes (in layers 4 and 5). Non-neuronal cells were largely seen in layer 1 and white matter.

To discover yet-finer volumetric patterns, we further analyzed the distribution of distances from each individual cell from each sequencing-defined subtype to its nearest neighbors, finding unexpectedly that the nearest neighbor of any inhibitory neuron tended to be its own subtype, rather than excitatory neurons or other inhibitory subtypes (Fig. 6G). If inhibitory neurons were randomly dispersed among the more abundant excitatory neurons in a purely salt-and-pepper distribution, the distance between inhibitory neurons would be larger than that from inhibitory to excitatory neurons (Fig. 6H). Instead, the actual intrasubtype distance of inhibitory neurons was much shorter (~15 μm, which is equivalent to the size of a single neuron, indicating direct somatic juxtaposition) (Fig. 6I), revealing a short-range self-clustering organization of inhibitory subtypes across volumes that could only be accurately measured in three but not in two dimensions (fig. S18A). When guided by this initial STARmap observation, evidence for such patterning could be also obtained in transgenic mouse lines (fig. S18, B and C). This discovery bears considerable relevance to previous functional work; for example, electrophysiological studies have suggested that inhibitory neurons in spatial proximity tend to be connected by electric (gap) junctions important for setting up synchronized firing patterns (40, 41), and in vivo imaging has suggested that inhibitory-neuron groupings in visual cortex could sharpen visual responses (42).

Discussion

STARmap defines a platform for 3D in situ transcriptomics, enabled by state-of-the-art DNA library preparation and sequencing in an HTC formulation. Here, STARmap was shown to be applicable to the study of molecularly defined cell types and activity-regulated gene expression in mouse cortex and to be scalable to larger 3D tissue blocks so as to visualize short- and long-range spatial organization of cortical neurons on a volumetric scale not previously accessible. In future work, STARmap may also be adapted to longer sequencing lengths or higher gene numbers; there is no intrinsic limit to the number of genes or RNA species that can be simultaneously and quantitatively accessed with STARmap (fig. S15); STARmap may also be capable of integrating cell-type information with single-neuron morphology and projection anatomy (for example, by means of Brainbow and MAPseq) (43, 44) as well as with in vivo neural activity imaging and electrophysiology. This platform can also be generalized to study other heterogeneous cell populations in diverse tissues across the body, although the brain poses special challenges well suited to STARmap analysis. For example, the polymorphic ARG expression observed across different cell types is likely to depend on both intrinsic cell-biological properties (such as signal transduction pathway-component expression) and on extrinsic properties such as neural circuit anatomy that routes external sensory information to different cells (here, in visual cortex). In general, it may not be possible to fully define brain cell typology independent of such 3D anatomy as well as activity patterns exhibited and experienced by cells during behavior; the nature of input and output communication pathways for the cells in question in fact can form the foundation for defining cell types (35, 36). Toward this end, in situ transcriptomics exemplified by STARmap can effectively link this imaging-based molecular information with complementary cellular-resolution datastreams describing anatomy, natural activity, and causal importance, thus promising to fundamentally deepen our understanding of brain function and dysfunction (2).

Methods summary

All animal procedures followed animal care guidelines approved by Stanford University’s Administrative Panel on Laboratory Animal Care (APLAC) and guidelines of the National Institutes of Health. For thin sections, animals were anesthetized and rapidly decapitated; the brain tissues were sliced by use of a cryostat. For thick sections, animals were anesthetized and transcardially perfused with paraformaldehyde; the brain tissues were sliced by use of a vibratome. In STARmap experiments, tissues were hybridized with SNAIL probes, enzymatically amplified, hydrogel embedded, and sequentially imaged by using the SEDAL process and a confocal microscope. The resulting image datasets were registered across multiple cycles by using the positions of all amplicons in each cycle and decoded. For cell typing and single-cell gene expression analyses, the amplicons were attributed to individual cells based on segmentation images of fluorescent Nissl staining. All the detailed procedures for the experiments and data analyses are described in the supplementary materials.

Supplementary Materials

www.sciencemag.org/content/361/6400/eaat5691/suppl/DC1

Materials and Methods

Figs. S1 to S18

References (4549)

Tables S1 and S2

References and Notes

Acknowledgments: We are grateful for the efforts and support of Deisseroth laboratory members A. Crow, S. Quirin, A. Chibukhchyan, C. Lee, M. Lo, and C. Delacruz. We thank M. Lovett-Barron for helping with image registration and J. H. Lee (Cold Spring Harbor Laboratory), J. Yan (Princeton University), N. D. Donoghue (Brown University), and Q. Dai (University of Chicago) for advice. We thank Deisseroth laboratory members L. Ye, A. Andalman, H. Wang, and B. Hsueh for discussion. We also thank L. Luo for suggestions on the manuscript. Funding: X.W. is supported by a Life Science Research Foundation fellowship and the Gordon and Betty Moore Foundation. W.E.A. is supported by a Fannie and John Hertz Foundation Fellowship and an NSF Graduate Research Fellowship. E.L.S. is supported by a NIMH Ruth L. Kirschstein NRSA fellowship (1F32MH110144-01). M.A.W. is supported by an NIMH Career Development Award (1K08MH113039). J.L. is supported by a Bio-X Interdisciplinary Initiatives Seed Grant. G.P.N., N.S., and F.-A. B. were supported by the Parker Institute for Cancer Immunotherapy, the FDA, and the NIH; F.-A.B. was supported by the Human Frontiers Science Program. K.D. is supported by NIMH (R01MH099647), NIDA (P50DA042012), the DARPA NeuroFAST program W911NF-14-2-0013, the NSF NeuroNex program, the Gatsby Foundation, the AE Foundation, the NOMIS Foundation, the Fresenius Foundation, the Wiegers Family Fund, the James Grosfeld Foundation, the Sam and Betsy Reeves Foundation, and the H. L. Snyder Foundation. Author contributions: X.W. and K.D. initiated the STARmap project to integrate HTC with in situ sequencing. X.W. developed the STARmap HTC, SEDAL sequencing, and STARmap hardware, designed all the DNA probes,and conducted the experiments. W.E.A. developed the STARmap software pipeline and analyzed the sequencing data. N.S., F.-A.B., and G.P.N. designed the initial versions of the SNAIL system. X.W. and F.-A.B. developed the SNAIL process for mouse brain tissue and compared SNAIL with other in situ methods. M.A.W. and E.L.S. conducted animal behavior and preparation of mouse brain tissue and contributed valuable advice. M.A.W. compared Nissl staining with other segmentation methods. W.E.A. and E.L.S. designed the visual-stimulus and cocaine-stimulus procedures. S.V. generated the CLARITY data with PV transgenicmice. K.E. and C.R. contributed to preparation of cell cultures. C.L. assisted with experiments. J.L. assisted with experiments and graphic design. K.D. supervised all aspects of the work. X.W., W.E.A., and K.D. interpreted the STARmap data, designed and prepared the figures, and wrote the manuscript, with edits from all authors. Competing interests: The design, steps, and applications of STARmap are covered in pending patent application material from Stanford University; all methods, protocols, and sequences are freely available to nonprofit institutions and investigators. Data and materials availability: All data are available in the main text or the supplementary materials. Computational tools and code and other materials are available at http://clarityresourcecenter.org. References for the open-source software packages used here can be found in the supplementary materials.
View Abstract

Navigate This Article