Sensory coding mechanisms revealed by optical tagging of physiologically defined neuronal types

See allHide authors and affiliations

Science  13 Dec 2019:
Vol. 366, Issue 6471, pp. 1384-1389
DOI: 10.1126/science.aax8055

Cell-type identification in neural circuits

In many cases, a single molecular marker is insufficient to define a specific cell type and may label a few, or a few hundred, physiologically distinguishable cell types. Lee et al. developed a high-throughput technique, called physiological optical tagging sequencing (PhOTseq), for identifying the expression profile of cells that exhibit a particular physiological profile (see the Perspective by Renninger). They used PhOTseq to identify genes encoding vomeronasal receptors in mice, which detect pheromones and subserve social communication.

Science, this issue p. 1384; see also p. 1311


Neural circuit analysis relies on having molecular markers for specific cell types. However, for a cell type identified only by its circuit function, the process of identifying markers remains laborious. We developed physiological optical tagging sequencing (PhOTseq), a technique for tagging and expression profiling of cells on the basis of their functional properties. PhOTseq was capable of selecting rare cell types and enriching them by nearly 100-fold. We applied PhOTseq to the challenge of mapping receptor-ligand pairings among pheromone-sensing neurons in mice. Together with in vivo ectopic expression of vomeronasal chemoreceptors, PhOTseq identified the complete combinatorial receptor code for a specific set of ligands.

Molecular markers are a powerful tool for labeling and analyzing neuronal cell types. However, in many cases, a single marker is insufficient to define a distinct cell type and may label a few, or a few hundred, physiologically distinguishable cell types (13). In such cases, one wishes to select specific physiological populations and discover their molecular identities. However, tools for proceeding from function to molecular markers are not fully mature. Markers such as c-fos label active neurons promiscuously across many cell types (4). Cells can be chosen for sequencing based on patch recording (58); however, this approach faces obstacles for rare cell types or when one needs many cells of the same type to profile low-abundance transcripts.

We reasoned that the key bottlenecks of this pipeline could be replaced by an all-optical approach. One prototype is CaMPARI, which uses the simultaneous presence of calcium and excitation by short-wavelength light under the control of the experimenter (9). Although CaMPARI allows one to control the labeling of excited neurons temporally, as with c-fos, any active neuron will be labeled, so this may comprise dozens or more cell types. There is considerable need for tools that provide the specificity of patch-seq with the convenience and high throughput of optical methods.

Here, we report a method called physiological optical tagging sequencing (PhOTseq). PhOTseq allows specific cell types to be defined intersectionally through their physiological activity under diverse conditions. We exploited two coexpressed fluorescent proteins, GCaMP and photoactivatable mCherry (PAmCherry) (1012); GCaMP permits recording neuronal activity by large-scale calcium imaging, and PAmCherry enables long-term labeling of selected neurons by photoactivation (PA). Using these reporters, we performed calcium imaging with a diverse set of stimuli or experiences and then selected a specific collection of cells for tagging by PA. Last, the tagged neurons were harvested and profiled for their mRNA expression (Fig. 1A).

Fig. 1 Two-photon PA can tag neurons chosen by activity pattern.

(A) Workflow of PhOTseq. NA, numerical aperture; qPCR, quantitative polymerase chain reaction; RNASeq, RNA-sequencing; Stim, stimulus. (B to D) Two-photon calcium imaging of the VNE explanted from a GCaMP5g-p2a-PAmCherry transgenic mouse. (B) Representative images of calcium chemosensory response evoked by either 10 μM A7864 (top left) or 10 μM E1050 (top right). PA mask (bottom) selects cells responsive to both stimuli. Scale bar, 20 μm. (C) The average GCaMP intensity obtained from masked cells. Black bars: delivery time course of stimuli. (D) Individual responses (ΔF/F) among masked cells to 10 μM A7864 and 10 μM E1050. (E) Confocal imaging after online two-photon PA. PAmCherry signals (left), PAmCherry signals with two-photon PA mask (middle), and PAmCherry signals with calcium responses to either 10 μM A7864 (cyan) or 10 μM E1050 (green) (right). Red and yellow arrow heads: the examples of PA mask regions associated with single and double PAmCherry-positive neurons, respectively. Scale bar, 20 μm. (F and G) FACS analysis of cells from nonphotoactivated tissue (F) or photoactivated tissue (G). Clusters were marked as P1, P2, P3, P4, and P5. P1 and P2 represent a photoactivated and an experimental control group, respectively. Of the total population, P1 and P2 accounted for 0.013% and 33.9% (F) or 0.218% and 37.3% (G), respectively. AU, arbitrary unit.

We created several plasmids and viruses, all expressing variants of GCaMP-2A-PAmCherry. The self-cleaving 2A peptide (13) allowed both proteins to be expressed from a single mRNA; the tight coupling facilitated ratiometric analysis to compensate for variability in expression. Initial in vitro and in vivo experiments demonstrated the feasibility of performing both calcium imaging and phototagging by using this reporter (figs. S1 and S2). To test the performance of the complete pipeline (Fig. 1A), we turned to mouse vomeronasal sensory neurons (VSNs), whose primary function is pheromone sensing and social communication. VSNs exhibit dense cell packing (14)⁠ (thus challenging phototagging accuracy), extreme functional heterogeneity, and an unambiguous relationship between molecular identity and physiological function (1518). We generated a transgenic tetO-GCaMP5g-2A-PAmCherry mouse and drove expression in all VSNs through olfactory marker protein (OMP)–internal ribosomal entry site (IRES)–tetracyclin transactivator (tTA) (fig. S3A). We observed calcium responses from the VSNs, and two-photon PA resulted in phototagging at single-cell resolution (figs. S3 and S4).

We evoked combinatorial VSN responses using two ligands: 5-androsten-3β, 17β-diol disulfate (A7864) and 1,3,5(10)-estratrien-3, 17β-diol disulfate (E1050). A custom online algorithm automatically segmented responsive neurons; out of 430 neurons responding to at least one ligand, 192 responded to both ligands. Voxels corresponding to these dually responsive neurons were chosen to create a mask so that only these voxels were illuminated during PA (Fig. 1, B to D). Photoactivated neurons were readily distinguished from the background, and PA tagged an average of 1.32 cells per region (Fig. 1E and fig. S5). The photoactivated vomeronasal epithelia (VNEs) were subsequently subjected to dissociation and fluorescence-activated cell sorting (FACS) (Fig. 1, F and G). In the FACS analysis, four to five abundant clusters were observed, of which only one (P1) was almost entirely dependent on PA. This population represents the one selected later for sequencing.

We used PhOTseq to begin unraveling the logic of chemosensation. Fewer than a dozen vomeronasal receptors (VRs) have ligands of known structure (5, 19, 20), and these are scattered widely across the gene family. We reasoned that a saturation analysis of VRs responsible for encoding a focused region of “chemical space” will provide insights into the molecular logic of chemosensation. First, we examined functional types of VSNs broadly by light-sheet calcium imaging while presenting 15 different sulfated steroids, of which six were also delivered over a range of concentrations, for a total of 27 different chemosensory stimuli. Approximately 15,000 neurons, ∼10% of the total VSN population per hemisphere, were visualized in each imaging volume (3). From three imaging volumes, more than 5000 neurons responded to at least one of the ligands, and they fell naturally into 20 different response types (Fig. 2, A and B). From this reference clustering, we chose to phototag the most abundant cell type, cluster 1, and three others showing similar chemoreceptive fields (clusters 2, 3, and 6) (fig. S6), among which cluster 3 comprised fewer than 2% of the responsive neurons and presumably ∼0.2% of the total VSN population (Fig. 2B). As a control group, we selected one of the most dissimilar cell types, cluster 19, which strongly responded to sulfated pregnanolones.

Fig. 2 PhOTseq-identified VR genes of overlapping cell types.

(A) Example calcium traces (normalized ΔF) from two cells. Amplitude was normalized by the maximum amplitude of all cells recorded. Ligands are listed at the top. OCPI, objective coupled planar illumination. (B) Neuronal responses to sulfated steroids. Cells are on columns, and stimuli are on rows. If not indicated, the ligand concentration is 10 μM. The color bar indicates normalized response. Cluster identities are reported at the top. Asterisks indicate PhOTseq target cell types. Asterisks in parentheses indicate functional type whose receptor identity was discovered during analysis by ectopic expression in Fig. 3. (C) Expression of the 30 most highly expressed VR genes when averaged across all sequenced cells; also shown are three marker genes (Omp, Gnai2, and Gnao1). Cells are on columns, and genes are on rows. PhOTseq-targeted functional types are shown at the top; cells belonging to these functional clusters were specifically photoactivated and sequenced. The color bar indicates the log-normalized expression level. ID, identification. (D) Proportion of a VSN type in each group shown in (C). Each tick on the x axis represents a different VR gene. Each functional type exhibited only one common VSN type. (E) Expression of five VR genes across different experiment groups. “P” indicates photoactivated cells, and “C” indicates nonphotoactivated control cells. *Padj < 0.01 (Wilcoxon rank sum test; Padj < 0.01) (Wilcoxon rank sum test; Padj = 2.7 × 10−11, 2.4 × 10−10, 1.1 × 10−20, 1.3 × 10−40, and 2.7 × 10−30 ; average fold difference: 6.2, 23.8, 80.7, 75.7, and 45.9 for Vmn1r89, Vmn1r86, Vmn1r78, Vmn1r237, and Vmn1r58, respectively). (F) Single-cell expression of Vmn1r85 and Vmn1r86 is nonexclusive.

For cell type selection under two-photon microscopy, on the basis of Fig. 2B, we used a subset of stimuli that distinguish cells among these five selected clusters (figs. S7 and S8 and movie S1). After online two-photon PA, we collected photoactivated or nonphotoactivated experimental control cells by FACS and performed single-cell RNA sequencing. We obtained a total of 622 qualified cells (fig. S9), among which 61, 44, 49, 40, and 65 photoactivated cells were from experiments aimed at functionally defined cluster 1, cluster 2, cluster 3, cluster 6, and cluster 19 cells, respectively (Fig. 2C).

Physiological responses of chemosensory neurons are thought to be largely determined by the VR genes they express (21). Therefore, our investigation was focused on the expression of VR genes. Among photoactivated cells, any given cluster exhibited substantial enrichment of at least one VR gene, with different VR genes being enriched in different clusters (Fig. 2C). By contrast, among control cells, VR gene expression was sporadic (fig. S10). Because our targeting accuracy was less than 100% (fig. S5C), we also expected some sporadic expression among cluster-selected photoactivated neurons (Fig. 2C). When the type of a VSN was defined by its maximally expressed VR gene, only one abundant type was found among cluster-selected photoactivated neurons (Fig. 2D) [cluster 1: Vmn1r89 type (25%); cluster 2: Vmn1r86 type (30%); cluster 3: Vmn1r78 type (27%); cluster 6: Vmn1r237 type (50%); cluster 19: Vmn1r58 type (58%)]. These five VR genes were the only significantly enriched receptor gene in each photoactivated group compared with all the other cells (Fig. 2E). During these analyses, we updated the Vmn1r237 gene model because the read coverage of Vmn1r237 and cloning revealed 3′ untranslated regions missing from existing gene models (fig. S11 and data file S1). We also unexpectedly observed coexpression of Vmn1r85 in Vmn1r86 neurons (Fig. 2F and fig. S12, A to C). Including data from our control and sporadic cells, coexpression was also observed for several other VR gene pairs, in each case consisting of a genomically adjacent pair (fig. S12, D to F). Taken together, we identified six putative VR genes mediating a focused set of responses.

To test whether PhOTseq identified true receptor-ligand pairs, we performed a gain-of-function study. Because of abnormal localization of VR proteins (22), in vitro heterologous expression systems have had little success, particularly for the V1R family (23). We reasoned that a VSN cell’s endogenous machinery would allow functional expression of a VR gene; seeking a more efficient route than making a transgenic mouse (5), we explored a virus-mediated approach. Among the viral vectors we tested, intravenous injection of rAAV2/8-CAG–green fluorescent protein (GFP) was able to induce expression of GFP in a subset of VSNs (fig. S13). We delivered rAAV2/8-CAG-GCaMP-2A-VR, in which VR was one of the genes identified by PhOTseq, and for each obtained multiple neurons that showed statistically significant responses to at least one ligand (Fig. 3, A and B, and movie S2). In each case, the dominant response pattern matched the expected PhOTseq response pattern (Fig. 3, C and D). We also tested Vmn1r85 because this gene was found to be coexpressed with Vmn1r86 (Fig. 2C). The response pattern matched that of cluster 4, which was another group that was similar to cluster 1 but was not chosen for analysis by PhOTseq (Fig. 3, C and D). Because sequencing revealed control neurons that expressed only Vmn1r85 (Fig. 2F), functional cluster 4 cells likely expressed Vmn1r85.

Fig. 3 Ectopic expression–enabled functional analysis of VRs.

(A) Expression by means of AAV injection in the temporal vein of newborn pups. AAV, adeno associated virus; P0.5, postnatal day 0.5. (B) Optical section from light-sheet calcium imaging after the ectopic expression of GCaMP-2A-Vmn1r237 in response to different ligands. The apical layer is occupied by dendritic tips (Den) with cell bodies (CB) below. Scale bar, 20 μm. (C) Calcium response after ectopic expression of GCaMP-2A-Vmn1r89, -Vmn1r86, -Vmn1r78, -Vmn1r237, -Vmn1r58, or -Vmn1r85. Ligands are identical to those in Fig. 2B. (D) Pairwise correlation between the reference clustering and the ectopic responses shown in (C) (left) or autocorrelogram of the reference clustering (right).

Only a few studies have examined the relationship between VR sequence and chemosensory function (5, 19) and never from the vantage point of complete knowledge of the VRs that underlie a set of nearest neighbors in terms of function. To investigate this relationship, we performed multiple sequence alignment (MSA) of vomeronasal type 1 receptor (V1R) protein sequences followed by phylogenetic tree analysis (Fig. 4A). Among V1Rs expressed in functionally similar types, VMN1R89, VMN1R86, and VMN1R85 were close to each other in the sense of putative evolutionary distance. Two of the functionally similar cell types, VMN1R78 and VMN1R237, belonged to different branches. In this representation, VMN1R58, which was the chosen outlier in terms of function, was not notably more divergent from VMN1R89, VMN1R86, and VMN1R85 than the more functionally similar VMN1R78 and VMN1R237. Consequently, this tree structure representation did not suggest a particularly close correspondence between function and primary sequence.

Fig. 4 Sequence and chemoreceptive similarity are strongly correlated.

(A) Unrooted phylogenetic tree of V1R proteins. The five functionally similar (blue) and one distant (red) receptor studied here are marked. The scale bar indicates the number of amino acid substitutions per site. Numbers indicate bootstrap values. (B) Pairwise distances among 204 V1R protein sequences (see supplementary materials, materials and methods). (C) For these deorphanized VRs, the top 20 nearest VRs, based on (B), are rank ordered. The deorphanized VRs are colored if shown in the table. (D) Two-dimensional representation of the distance matrix shown in (B). Each dot represents a single VR protein. coord, coordinate.

However, when the pairwise distances acquired from the MSA analysis were examined carefully, we noticed that, despite their apparent evolutionary distance, VMN1R78 and VMN1R237 were consistently among the VR genes most closely related to VMN1R89, VMN1R86, and VMN1R85 (Fig. 4, B and C). Among 204 V1R sequences examined, the median neighbor rank among these VRs was 11, implying that these VRs are typically among the top-few-percent–closest matches to one another. We therefore considered the possibility that phylogeny, in attempting to construct plausible ancestry, inaccurately models relatedness between sequences that are not nearest neighbors. To more comprehensively evaluate and visualize all pairwise distances, we performed a classical multidimensional scaling analysis to project the high-dimensional sequence-based distance relationships into two dimensions (Fig. 4D). This resulted in a different picture of the gene family, whose dominant feature was the presence of seven to nine apparent clusters. The functionally similar receptors, including VMN1R78 and VMN1R237, were near neighbors, whereas the functionally dissimilar receptor VMN1R58 was positioned in a distant cluster. We conclude that relationships among the primary sequences of these VRs are strongly correlated with their degree of functional similarity.

Our results demonstrate the utility of PhOTseq, an all-optical solution to the problem of cell type selection and labeling. It provides an opportunity to study the relationships between genes and physiological function even in extremely rare cell types that can be defined only through extensive functional characterization. PhOTseq will be applicable when individual neurons respond combinatorially to various conditions, including sensory stimulation, emotional status, and behavior. We further anticipate that in vivo ectopic expression will be widely applicable for characterizing VR-ligand pairings. In addition, our data show unexpected receptor coexpression, a specific and systematic departure from the “one neuron, one receptor rule” (21)⁠, whose potential roles will need to be evaluated in many aspects of olfactory function. Lastly, our study provides biological insight on chemosensation by comprehensively mapping receptor-ligand pairings for chosen subsets of the vomeronasal sensory population and suggests that saturation analyses can reveal sequence-function coupling unanticipated by single-ligand studies (24, 25).

Supplementary Materials

Materials and Methods

Figs. S1 to S13

References (2652)

Movies S1 and S2

Data File S1

References and Notes

Acknowledgments: We thank R. Yu for sharing the TetO vector and OMP-IRES-tTA mouse, P. H. Taghert for the HEK cell line, H. Schoknecht for mouse care, and D. W. Kim for developing imaging software. We also thank T. Barnes, J. Chen, J. Corbo, X. Fu, C. Greer, A. Livi, R. Roberts, and M. Zhang for suggestions and comments. This work was supported by the Hope Center Viral Vectors Core and the Mouse Genetics Core at the Washington University School of Medicine in St. Louis. Funding: This work was funded by NIH/NIDCD R01 DC005964 and DC010381. Author contributions: M.K. performed in vivo imaging and PA. D.L. performed all other experiments and analysis. T.E.H. supervised the study and analysis. Competing interests: The authors have no competing interests. Data and materials availability: Accession numbers are available in the supplementary materials.

Stay Connected to Science

Navigate This Article