Physiology-forward identification of bile acid–sensitive vomeronasal receptors

See allHide authors and affiliations

Science Advances  29 May 2020:
Vol. 6, no. 22, eaaz6868
DOI: 10.1126/sciadv.aaz6868


The mouse accessory olfactory system (AOS) supports social and reproductive behavior through the sensation of environmental chemosignals. A growing number of excreted steroids have been shown to be potent AOS cues, including bile acids (BAs) found in feces. As is still the case with most AOS ligands, the specific receptors used by vomeronasal sensory neurons (VSNs) to detect BAs remain unknown. To identify VSN BA receptors, we first performed a deep analysis of VSN BA tuning using volumetric GCaMP6f/s Ca2+ imaging. These experiments revealed multiple populations of BA-receptive VSNs with submicromolar sensitivities. We then developed a new physiology-forward approach for identifying AOS ligand-receptor interactions, which we call Fluorescence Live Imaging for Cell Capture and RNA sequencing, or FLICCR-seq. FLICCR-seq analysis revealed five specific V1R family receptors enriched in BA-sensitive VSNs. These studies introduce a powerful new approach for ligand-receptor matching and reveal biological mechanisms underlying mammalian BA chemosensation.


Chemosensory systems extract salient information from environmental cues that are critical for social and reproductive behaviors. In mice and other terrestrial mammals, the accessory olfactory system (AOS) guides innate behaviors such as territorial aggression and mating (1, 2). The initial detection of olfactory stimuli in the AOS is mediated by vomeronasal sensory neurons (VSNs), located in the vomeronasal organ (VNO), before these signals are routed to the accessory olfactory bulb and then behaviorally important regions including the medial amygdala and bed nucleus of the stria terminalis.

Although the behavioral impacts of AOS chemosensation are manifold, our knowledge of the full complement of natural ligands for VSNs remains incomplete. Natural AOS ligands are found in the excretions of conspecific and heterospecific animals (e.g., urine, feces, tears, and saliva). These natural ligands include, but are not limited to, polar steroids, bile acids (BAs), major urinary proteins, and major histocompatibility complex peptide ligands (36). As we learn more about the full natural complement of AOS ligands, new questions are arising about how these cues are detected by specific types of VSNs. This knowledge is critical as we seek to understand how patterns of VSN activation by complex blends of environmental chemosignals cause changes in animal physiology and behavior.

A major barrier preventing a deeper understanding of AOS function is a lack of ligand-receptor matches for AOS cues. VSNs generally express just one or two of ~300 unique vomeronasal receptors (VRs) or formyl peptide receptors (1, 7). The largely monoallelic expression of individual VRs by VSNs means that each VSN takes on the chemosensory sensitivity (or “tuning”) of the dominantly expressed receptor. In theory, this feature would greatly simplify the process of unambiguously matching VSN ligands to their cognate receptors. However, efforts to generate large-scale VR screening assays via heterologous expression has, for the most part, been unsuccessful [although see (8)]. Other efforts have explored VNO ligand-receptor interactions using various approaches (912). Despite recent progress, a small fraction of the VR family has confirmed ligands (10, 11, 1316).

Given the importance of AOS ligand-receptor interactions for mammalian physiology and social behavior, we sought to develop a pipeline for matching AOS ligands to receptors. We devised a physiology-forward approach using native VNO tissue (with endogenous VR expression) to identify the receptors for four recently discovered natural BA ligands present in feces (4). BAs are steroids that are structurally similar to other potent AOS ligands including sulfated and carboxylated glucocorticoids (6, 17). Despite evidence that many VSNs show selectivity for BAs compared to these other steroids and evidence that some VRs are exquisitely sensitive to other steroids (911, 17), questions remain about whether there are BA specialist VRs or whether BAs nonselectively activate certain steroid-sensitive VRs.

Here, we report the discovery of distinct populations of VSNs that are BA specialists and others that are BA generalists. We comprehensively mapped the sensory space of VSNs to four prominent BAs and several sulfated steroids, including sex steroids and glucocorticoids, finding that some BA generalists are also sensitive to sulfated glucocorticoids. We report the successful development of a function-forward strategy, which we term Fluorescence Live Imaging and Cell Capture for RNA-seq (FLICCR-seq), to isolate highly BA-sensitive VSNs for single-cell RNA sequencing. We used FLICCR-seq to show that BA-sensitive VSNs are enriched in the expression of at least five V1R family VRs, two of which appear to be associated with broad BA tuning, and three that are associated with selective detection of specific BAs. Collectively, these studies reveal mechanisms of mammalian BA chemosensation and improve our understanding of chemosensory information as it enters the behaviorally important AOS pathway.


Sensitive and selective VSN responses to BAs

BAs are naturally occurring AOS ligands that are found abundantly in mouse and other animal feces (4, 18). However, the sensitivity and selectivity of VSNs to BAs has not been deeply explored. To quantify VSN BA tuning, we performed population VSN Ca2+ imaging via objective-coupled planar illumination (OCPI) microscopy (Fig. 1A). We began by stimulating VSNs with a concentration series of two closely related monomolecular BAs at concentrations spanning 100 nM to 10 μM (Fig. 1B). Monomolecular BA ligands elicited reliable, concentration-dependent responses from VSNs across multiple randomized, interleaved trials in this preparation (Fig. 1, C and D). We used both GCaMP6s- and GCaMP6f-expressing preparations (VNO epithelia from both OMPxAi96 and OMPxAi148 mice) throughout this manuscript.

Fig. 1 Evaluating VSN responses to monomolecular BAs with population VSN Ca2+ imaging via OCPI microscopy.

(A) OCPI (light sheet) microscopy imaging setup. Thousands of VSNs in the intact vomeronasal epithelium were imaged via this setup to measure the responses across a panel of ligands. (B) Experimental stimulus panel consisting of five test concentrations of two BAs (CA, cholic acid; DCA, deoxycholic acid) along with positive controls (FF, dilute female mouse feces) and negative controls (Ringer’s). (C) Representative fluorescence intensity traces of volumetric regions of interest (ROIs) isolating two individual VSNs over the course of one experiment. Colored bars represent ligand presentation corresponding to ligands depicted in (B). Each stimulus presentation lasted five stacks (~15 s), with 15 stacks (~45 s) between stimulus presentations. (D) Representative across-trial responses of the same VSNs in (C). Bolded traces show across-trial mean waveforms, with red line color indicating a statistically significant increase compared to Ringer’s control trials.

From these datasets, we first performed a pairwise comparison between VSN responses to a primary BA, cholic acid (CA), and its gut bacteria-dependent derivative deoxycholic acid (DCA). These molecules differ from each other only by the presence or absence of a hydroxyl group on the 7-carbon (Fig. 2). Previous studies of VSN activity using OCPI microscopy estimated that, for a larger panel of ligands, stimulus-responsive VSNs represented 20 to 50% of all VSNs (19). We use a more focused panel of ligands here and analyzed only well-registered VSNs showing spontaneous or stimulus-driven activity. We estimate that the analyzed fraction of VSNs in these studies represents ~25% of all VSNs. Within this active fraction of analyzed cells, we observed many VSNs that reliably responded to one or both of the BA ligands at 10 μM or lower (Fig. 2, A and B). We observed VSN populations that responded selectively to CA only (16.01%, corresponding to ~4% of all VSNs), DCA only (16.58%, ~4% of all VSNs), and a population that responded to both CA and DCA (20.21%, ~5% of all VSNs; Fig. 2C). As a control for possible nonselective effects of BAs on VSN Ca2+ signaling and effects of VSN spontaneous activity (movie S1), in this and all experiments, we analyzed responses of VSNs that demonstrated a visible change in GCaMP6 fluorescence during at least one stimulus trial. We found that cells chosen based on single-trial responses (i.e., including spontaneously active cells) were not reliably responsive to either CA or DCA (47.2%), suggesting that BAs do not cause nonselective VSN Ca2+ increases. Within the subset of VSNs that responded selectively to CA, a small but consistent subset (1.2% of analyzed VSNs, corresponding to ~0.3% of all VSNs) was sensitive to CA at 0.1 to 0.3 μM (Fig. 2C), indicating that—for this comparison—these neurons were CA specialists. The ~0.3% (~1 of 300) value represents the fraction of VSNs that would express a single receptor if all of the ~300 VRs were expressed by equal fractions of VSNs. These experiments made it clear that VSNs, both individually and as a population, can distinguish between functionally and structurally similar BAs.

Fig. 2 VSN responses to CA and DCA.

(A and B) Representative colorized images of VSNs response (ΔF/F) in a single frame of and OCPI image stack (left). The responses to 10 μM CA (red), 1 μM CA (green), and 0.1 μM CA (blue) are shown in (A). The responses of VSNs to 10 μM CA (red), 1 μM CA (green), and 10 μM DCA (blue) are shown in (B). Across-trial VSN responses are plotted as individual traces. Bolded trace indicates the mean response across all stimulus repeats. Responses in the cyan box correspond to the VSN indicated by the cyan arrow, those in the orange box correspond to those indicated by orange arrows, and traces in the magenta box correspond to the neuron indicated by the magenta arrow. (C) Clustered heat map of VSN response to CA and DCA at varying concentrations. Each column indicates an individual neuronal response. Cluster divisions are indicated by black vertical lines. Arrows above the heat map highlight the clusters in which neurons from (A) and (B) fall. Shown are the responses of 1222 VSNs from three OMPxAi96 animals. All experiments included at least three randomized, interleaved stimulus repeats. (D) Chemical structure diagrams of CA and DCA. Red annotations indicate key differences. (E) Heat map of cluster separation. The map indicates the pairwise discriminability index (d′) between each cluster [based on linear discriminant analysis (LDA)]. White hue indicates the P = 0.05 cutoff, and shades of red indicate higher discriminability (lower P values). For all these pairwise comparisons, P < 0.05.

CA and DCA are both found in mouse feces and are closely related in structure, but previous work indicated that mouse VSNs are also sensitive to lithocholic acid (LCA), a BA that is derived from a different primary BA [chenodeoxycholic acid (CDCA)] and which is not present at detectable levels in mouse feces (Fig. 3) (4, 18). In separate experiments, we examined the response profiles of VSNs to CA and LCA across the same 100 nM to 10 μM concentration range (Fig. 3). Confirming earlier studies in the downstream AOB (4), we observed reliable VSN responses to LCA (Fig. 3, A and B). Clustering based off VSN response profiles revealed populations of neurons that responded selectively to CA only (21.5% of analyzed VSNs, ~5.4% of all VSNs), LCA only (14.8% of analyzed VSNs, ~3.7% of all VSNs), and both CA and LCA (11.0% of analyzed VSNs, ~2.75% of all VSNs; Fig. 3C). As before, a subset of CA-selective VSNs showed sensitivity to CA in the 0.1 to 0.3 μM range (0.45% of analyzed VSNs; Fig. 3C). In this experiment, we noted that many VSNs that were analyzed on the basis of spurious single-trial (likely spontaneous) activity showed substantial across-trial mean ΔF/F responses to 10 μM LCA (Fig. 3C). LCA can be cytotoxic (20), so the increased mean ΔF/F to 10 μM LCA in some VSNs may reflect mild cytotoxicity. Additional pairwise concentration-response experiments produced similar results: In each, we found evidence for BA specialists and generalists with submicromolar-to-micromolar BA sensitivities (figs. S1 to S4). Of these other comparisons, we found that CA and CDCA, both primary BAs, showed the highest degree of coactivation (fig. S1). Cumulatively, these data indicate that, as we previously hypothesized, there are likely to be many BA-sensitive receptors expressed by VSNs.

Fig. 3 VSN responses to CA and LCA.

(A and B) Representative colorized images of VSNs (ΔF/F) in a single frame of an OCPI image stack (left). The responses to 10 μM CA (red), 1 μM CA (green), and 10 μM LCA (blue) are shown in (A). The responses of VSNs to 10 μM CA (red), 1 μM CA (green), and 0.1 μM CA (blue) are shown in (B). Across-trial VSN responses are plotted as individual traces. Bolded trace indicates the mean response across all stimulus repeats. Responses in the cyan box correspond to the VSN indicated by the cyan arrow, while those in the orange box correspond to those indicated by orange arrows. (C) Clustered heat map of VSN response to CA and LCA at varying concentrations. Each column indicates an individual neuronal response. Cluster divisions are indicated by black vertical lines. Arrows above the heat map highlight the clusters in which neurons from (A) and (B) fall. Shown are the responses of 1362 VSNs from three OMPxAi96 animals. All experiments included at least three randomized, interleaved stimulus repeats. (D) Chemical structure diagrams of CA and LCA. Red annotations indicate key differences. (E) Heat map of cluster separation. The map indicates the pairwise discriminability index (d′) between each cluster (based on LDA). White hue indicates the P = 0.05 cutoff, and shades of red indicate higher discriminability (lower P values). For all these pairwise comparisons, P < 0.05.

VSNs show diverse responses to BAs and sulfated steroids

The high degree of overlap in VSN sensitivities seen in pairwise BA tuning comparisons suggested that a broader assay would be helpful in further determining the breadth of VSN populations (and presumably the VRs) that are sensitive to these four common BAs. Moreover, questions remain about whether BA-sensitive VSNs are also activated by other steroid ligands (6, 13, 17, 21). We therefore conducted Ca2+ imaging experiments in which we exposed the VSNs to CA, DCA, CDCA, and LCA (each at 1 μM and 10 μM) and four sulfated steroids (each at 10 μM). This panel of sulfated steroids included an androgen (A6940), an estrogen (E0893), a pregnanolone (P3865), and a glucocorticoid (Q1570), each of which has been shown to activate specific subsets of VSNs (Fig. 4) (19, 22). As expected, VSNs showed reliable responses to this ligand panel (Fig. 4, A and B).

Fig. 4 VSN responses to a panel of monomolecular ligands and sulfated steroids.

(A and B) Representative colorized images of VSN responses (ΔF/F) in a single frame of an OCPI image stack (left). The responses to 10 μM CA (red), 10 μM DCA (green), and 10 μM Q1570 (blue) are shown in (A). The responses of VSNs to 10 μM CDCA (red), 10 μM A6940 (green), and 10 μM Q1570 (blue) are shown in (B). Across-trial VSN responses are plotted as individual traces. Bolded trace indicates the mean response across all stimulus repeats. Responses in the cyan box correspond to the VSN indicated by the cyan arrow, those in the orange box correspond to those indicated by orange arrows, and traces in the magenta box correspond to the neuron indicated by the magenta arrow. (C) Clustered heat map of VSN response to a panel of monomolecular ligands. Each ligand is shown at left to highlight structural differences, with red annotations indicating key differences among BAs. Each column indicates an individual neuronal response. Cluster divisions are indicated by black vertical lines. Arrows above the heat map highlight the clusters in which neurons from (A) and (B) fall. Shown are the responses of 890 VSNs from three OMPxAi96 animals. All experiments included at least three randomized, interleaved stimulus repeats. (D) Chemical structure diagrams of all ligands used in this experiment. Red annotations indicate key differences among BAs. (E) Heat map of cluster separation. The map indicates the pairwise discriminability index (d′) between each cluster (based on LDA). White hue indicates the P = 0.05 cutoff, and shades of red indicate higher discriminability (lower P values). For all these pairwise comparisons, P < 0.05.

Cluster analysis of VSN tuning to this broad BA and sulfated steroid panel revealed several interesting features (Fig. 4C). First, we found that a substantial population of VSNs (Fig. 4C, clusters 1 and 2) were sensitive to all four BAs in this panel at 10 μM, indicating that many of the VSNs that were coactivated in pairwise comparisons are broadly BA tuned. Within the broadly BA-tuned group of VSNs were cells that were also activated by 10 μM of the sulfated glucocorticoid Q1570 (Fig. 4C, cluster 1). This experiment provided more evidence of BA sensitive and selective VSN subpopulations (Fig. 4C, clusters 3 to 15) and confirmed many previously studied sulfated steroid–tuned VSN subpopulations (Fig. 4C, clusters 16 to 21) (11, 19, 21, 22). Also noteworthy were clusters 6, 7, 11, and 14, which showed some degree of coactivation by BAs and sulfated steroid ligands. Overall, we found that 29.1% of analyzed VSNs (~7.3% of all VSNs) were activated by at least one BA in this panel (at 10 μM), 32.3% of analyzed VSNs (~8.1% of all VSNs) were activated by at least one sulfated steroid, and 16.9% of analyzed VSNs (~4.2% of all VSNs) were activated by at least one BA and one sulfated steroid. The specific patterns of BA and sulfated steroid activity indicate, at a minimum, that there are dozens of VSN populations that encode BA identity. Moreover, the emergence of several BA- and sulfated steroid–sensitive populations indicates that VSNs use receptors (presumably VRs) that are sensitive to steroids with substantially different chemical features and which are derived from distinct biochemical pathways in the emitter.

Previous studies have determined that steroid-sensitive VSNs express members of the V1R subfamily, which happen to project axons selectively to the anterior subregion of the AOB (10, 11, 13, 22). On the basis of the superficial chemical similarities between BAs and other steroids, and the partial overlap seen between BA- and sulfated steroid–responsive VSNs, we hypothesized that V1Rs are likely serving as BA receptors. As a first pass at identifying the BA-sensitive VR family, we performed GCaMP6f Ca2+ imaging in specialized ex vivo preparations that maintain functional connections between the VNO and AOB (23). Using the same panel of BA and sulfated steroids ligands shown in Fig. 4, we stimulated the VNO in these ex vivo preparations while performing volumetric imaging of the AOB glomerular layer (where VSN axons terminate; fig. S5). Consistent with our hypothesis of V1R BA sensitivity, we observed glomerular layer activation in small, distributed regions of interest (ROIs) resembling individual glomeruli in the anterior AOB (fig. S5). Across these experiments, we did not identify a distinct spatial region or evidence of spatial colocalization between specific BA-tuned glomerular subpopulations (23), but more detailed investigations will be needed to fully evaluate the spatial distribution of steroid- and BA-sensitive glomeruli. Still, these data provide additional evidence that BA-sensitive VSNs express members of the V1R family. The accumulated evidence suggested that substantial populations of mouse VSNs are specifically tuned to naturally occurring BA ligands, implying that these VSN populations express sensitive and selective BA receptors.

Physiology-forward BA receptor identification via FLICCR-seq

The collective evidence of sensitive, selective, and diverse VSN BA tuning indicated that many BA-sensitive receptors are expressed by BA-sensitive VSNs. Identifying ligand-receptor interactions in the VNO remains a challenging task. We developed FLICCR-seq, a function-forward approach that uses GCaMP6 physiology at the front end, to identify and isolate BA-sensitive VSNs before single-cell RNA-seq (scRNAseq; Materials and Methods and Fig. 5A). We placed acute VNO slices from OMPxAi96 mice into the recording chamber of a modified slice electrophysiology rig. On the basis of OCPI Ca2+ imaging results, we selected 1 μM as the single BA concentration for stimulating BA-sensitive VSNs in these experiments (Figs. 2 to 4). As in OCPI experiments, we used computer-controlled stimulus delivery to apply no fewer than three randomized, interleaved ligand stimulations to each tissue to identify reliably BA-responsive VSNs (see Materials and Methods, Fig. 5, B and C, and movie S2). After we identified 1 μM BA-sensitive VSNs, we gently aspirated the stimulated cell into a modified patch pipette. The cell and its contents were then processed for scRNAseq. We selected a total of 158 samples, including 145 single VSN samples, seven bulk VNO samples, and six kit controls. Among the 145 single VSN samples were 27 CA-responsive, 16 DCA-responsive, 39 CDCA-responsive, and 22 LCA-responsive cells. We analyzed scRNAseq data using Seurat (24), finding gene expression patterns that included a larger, more dispersed group (clusters 0 to 2), a smaller, more compact cluster (cluster 3), and a tiny cluster representing kit (non-VNO) RNA controls (cluster 4; Fig. 5D). Inspecting the expression patterns of VSN-associated transcripts in these clusters revealed that the vast majority of samples in clusters 0 to 2 expressed Gnai2 (encoding the V1R-associated Gαi2 protein), whereas the vast majority of samples in cluster 3 expressed Gnao1 (encoding the V2R-associated Gαo protein; Fig. 5E). Other genes thought to be expressed by both V1R- and V2R-expressing VSNs, Trpc2 (encoding the TRPC2 cation channel) and Ano1 (encoding the anoctamin 1/TMEM16α calcium–activated chloride channel), were found to be expressed by the majority cells in all clusters except kit controls (cluster 4; Fig. 5E).

Fig. 5 Function-forward selection of BA-sensitive VSNs for scRNAseq.

(A) Overview of FLICCR-seq experimental setup. VNO slices from OMP-GCaMP6s mice were stimulated with 1 μM BAs a minimum of three times. BA-responsive cells are plucked with a modified glass patch pipette and processed for scRNAseq. (B) Representative fluorescence images of 1 μM CA-responsive VSNs. The VSNs marked with ‡ and † were selected and processed for scRNAseq. DIC, differential interference contrast. (C) Representative images of VSN selection during active stimulation (in this trial with 1 μM DCA). (D) UMAP multidimensional scaling of 158 samples, including 145 VSNs collected using these methods. Clusters 0 to 2 represent subsets of a larger group, whereas cluster 3 identifies a smaller, isolated cluster. Cluster 4 (gray, asterisk) contained RNA kit control samples and are accordingly distinct from VSNs. (E) Left: UMAP multidimensional scaling with each cell colorized based on expression of Gnai2 (Gαi, magenta) and Gnao1 (Gαo, cyan), genes that are associated with V1R- and V2R-expressing VRs, respectively. This analysis indicates that clusters 0 to 2 contain V1R-expressing VRs, while cluster 3 contains mostly V2R-expressing VRs. Right: Same UMAP scaling with each cell colorized based on the expression of Ano1 (Anoctamin 1, also known as TMEM16α) and Trpc2 (cyan), ion channels that have been studied in VSNs. Expression of Trpc2 or both Trpc2 and Ano1 was detected in almost all cells. (F) Genome viewer visualization of 18 representative VSNs showing the mapped scRNAseq reads for Trpc2, selected V1Rs and V2Rs, and G protein α subunits. The rows marked with ‡ and † reflect the cells indicated in (B) and (C). Note that the number of reads for individual VRs in each VSN was often stronger than for reference genes (e.g., Trpc2, Gapdh, and Actb).

Further inspection of the scRNAseq data for these samples revealed patterns of gene expression consistent with single VSNs, including the near exclusive expression of single VRs (Fig. 5F and fig. S6). In addition to the samples containing individual VSNs, we included a small number of VNO samples that included several hundred cells each. In these pooled samples, we detected the expression of many VRs at low expression levels (Fig. 5F and fig. S6), indicating that single-VR expression patterns were not an unintended feature of the scRNAseq pipeline, which necessarily involves several amplification steps. Collectively, the depth and quality of scRNAseq data were sufficient for further analysis to attempt to correlate gene expression patterns with patterns of BA responsivity for these VSN samples (Figs. 5 and 6).

Fig. 6 Identification of BA-sensitive VRs.

(A) Clustered heat map showing the gene expression for all 158 analyzed samples. Rows are individual transcripts, with divisions between functional subgroups indicated by a horizontal dotted line. At the top are control genes (see Materials and Methods). At the bottom is a lookup table of responsivity. Cyan indicates that the sample contained a VSN that was reliably responsive to the stimulus, dark blue indicates an inconsistent response, black indicates no response, and gray indicates that the ligand was not tested on that VSN. For cell amount (“Cell amt.”), green indicates a single cell was cleanly picked, orange indicates that additional cell material was visible in/on the pipette tip, red indicates a bulk tissue extraction, and black indicates RNA kit control samples. (B) Heat map plot of the top four differentially expressed genes for BA-responsive VSNs. In each subplot, the colored horizontal bars indicate the responsivity, with colors representing the same qualities as in (A). VSNs with unclear responsivity for the specified ligand were excluded from these plots. Expression levels (log counts/10,000) were normalized by the maximum values per gene across all samples. (C) Duplex RNA in situ hybridization micrographs for combined immediate-early gene (IEG) probes (Egr1 and Fos; blue stain) and test VRs (red stain). Mice were exposed in vivo to control saline or a mixture of BAs at 0.1 or 1 mM to stimulate VSNs. (D) Violin plot of the density of VR staining across VNO sections from 30 mice (15 males and 15 females). (E) Violin plot of the density of IEG staining (24 mice, three males and three females per condition, and four conditions). (F) Violin plots of the normalized number of colocalized IEG+ puncta per VR+ cell following in vivo exposure to BAs or female mouse fecal extracts. Multiple VNO sections were analyzed from 24 mice total [same cohort as in (E)]. †P < 0.1 and *P < 0.05 (Kruskal-Wallis test followed by multiple comparisons test). n.s., not significant.

Identification of five BA-sensitive V1Rs

Inspecting the expression patterns of individual scRNAseq samples revealed that several V1Rs were often expressed by BA-sensitive VSNs (Figs. 5F and 6A). For example, seven samples expressed Vmn1r16 (encoding the V1rc29 receptor), 14 samples expressed Vmn1r70 (encoding the V1rl1 receptor), and 36 samples expressed Vmn1r228 (encoding the V1re3 receptor; Figs. 5F and 6A). The enrichment of these VR genes in BA-sensitive VSNs was consistent with a correlation between the expression of these particular V1Rs and BA sensitivity. To determine the probability of achieving these results by chance, we performed differential gene expression analysis. First, we used responsivity to BAs as a grouping variable (Fig. 6B and table S1). Collapsing all BA-responses regardless of ligand identity (i.e., comparing all BA-responsive to all to BA-nonresponsive samples) revealed enrichment of Kirrel2 and Gnai2, which are known to be expressed by V1R-expressing VSNs (25, 26), but no strong correlation with any individual VR (Fig. 6B and table S1). When we included ligand identity as a grouping variable (e.g., comparing all 1 μM CA-responsive to all 1 μM CA-unresponsive samples), this analysis revealed enrichment of five individual V1Rs. Vmn1r16 (V1rc29) and Vmn1r228 (V1re3) were enriched in 1 μM CA-responsive samples, Vmn1r21 (V1rc28) was enriched in 1 μM DCA-responsive samples, Vmn1r225 (V1re5) was enriched in 1 μM CDCA-responsive samples, and Vmn1r70 (V1rl1) was enriched in 1 μM LCA-responsive samples (Fig. 6B and table S1).

We then used a separate approach for evaluating the statistical significance of these enriched receptors, this time by grouping samples based on their expression of the candidate V1Rs (Fig. 6A and table S2). Using this approach, we found that Vmn1r225- and Vmn1r228-expressing cells were statistically more likely to respond to any BA (regardless of BA identity) than expected by chance, indicating a potential role for these two V1re clade members in broader BA sensitivity (P < 0.01, binomial test; Fig. 6A and table S2). Vmn1r16-expressing cells were only enriched in 1 μM CA sensitivity, indicating that Vmn1r16 (V1rc29) may be a CA-selective receptor (Fig. 6A and table S2). Vmn1r21-expressing cells were enriched in 1 μM DCA sensitivity, indicating that Vmn1r21 (V1rc28) may be a DCA-selective receptor (Fig. 6A and table S2). Vmn1r70-expressing cells were enriched in 1 μM LCA sensitivity, indicating that Vmn1r70 (V1rl1) may be an LCA-selective receptor. Although many of the samples we collected expressed Vmn1r89 (V1rj2), which was previously shown to be highly sensitive to certain sulfated estrogens (13), we did not find evidence that Vmn1r89 or several other expressed receptors (e.g., Vmn2r1, Vmn2r2, and Vmn2r54) were enriched for BA responsiveness (table S2). Collectively, these analyses provide strong evidence for V1R enrichment in populations of VSNs with specific patterns of BA sensitivity and support the hypothesis that the diverse patterns of BA tuning observed via OCPI microscopy are due to the presence of several BA-sensitive VRs.

The statistical evidence for a correlation between expression of these V1Rs and BA sensitivity patterns was strong, so we sought additional verification of these ligand-receptor pairs. We exposed adult male and female mice in vivo to a mix of all four of the tested BAs at concentrations of 0.1 to 1 mM, which, given the expected 100-to-1000-fold dilution from the nares into the VNO lumen, would achieve effective BA concentrations of 0.1 to 10 μM. Following this in vivo stimulation, we performed RNA in situ hybridization for all candidate BA-sensitive VRs (Vmn1r16, Vmn1r21, Vmn1r70, Vmn1r225, and Vmn1r228), a negative control VR (Vmn2r54), and two pooled immediate-early genes (IEGs) (Egr1 and Fos), similar to previous studies (Fig. 6, C to F) (10). The number of VSNs expressing each VR was variable, with Vmn1r16 expressed by significantly fewer VSNs than all other VR probes (Fig. 6D). Labeling from the Egr1 + Fos probes increased with in vivo exposure to 1 mM BAs and female mouse feces extract (Fig. 6E). We observed increased colocalization between all candidate VR probes and our combined IEG probe (Egr1 + Fos) following in vivo BA stimulation at 0.1 to 1 mM (Fig. 6F). In contrast, we did not observe increased colocalization between IEG probes and the Vmn2r54 negative control (Fig. 6F). These experiments provided additional evidence that each of the candidate V1Rs identified by FLICCR-seq are sensitive to BA stimulation in vivo.


VSNs are sensitive, selective BA detectors

Previous studies showed that BAs can elicit robust AOS activation, and the diversity of BA tuning by individual neurons in the AOB (immediately downstream of VSNs) indicated that BAs may be sensed by multiple VSN receptors (4). BAs, which are found abundantly in animal feces, have several features that could provide useful information to the receiving animal. For example, BAs vary by species, sex, diet, and gut microbiome (18). Across the animal kingdom, BAs have long been studied as olfactory ligands and are used by several fish species as pheromones (27). Recent work in zebrafish has identified BA-sensitive ORA (olfactory receptor class A-related) class receptors in the olfactory epithelium (28). In the context of our growing knowledge of AOS ligands, BAs are a compelling but poorly understood source of mouse social chemosensory information. Hence, identifying the mechanisms of detection, projection, and processing of BA information will be important for producing a comprehensive view of the AOS and its many impacts on mammalian social and reproductive behaviors.

Although BAs were previously shown to reliably elicit AOS activation in the AOB (4) and VNO (29), there has been relatively limited information about the sensitivity and selectivity of VSNs to BAs. Other monomolecular AOS steroid ligands have been found to be highly potent VSN activators—some at subnanomolar concentrations—so we sought to thoroughly evaluate VSN BA tuning using OCPI microscopy (Figs. 1 to 4 and figs. S1 to S4). Volumetric Ca2+ imaging via OCPI microscopy remains one of the most robust methods for evaluating VSN ligand sensitivity because it enables the unbiased sampling of hundreds to thousands of VSNs per animal (19, 21, 29, 30). Concentration-response experiments using OCPI with pairs of BAs revealed the presence of multiple BA-sensitive VSN populations (Figs. 1 to 4 and figs. S1 to S4). Many of these VSNs showed submicromolar BA sensitivity, comparable to other steroid-sensitive VRs (9, 11).

These experiments show that each of the BAs tested is capable of selectively driving a subset of VSNs. However, some questions remained about whether VNO BA sensitivity was a by-product of sensitivity to other steroids (e.g., polar sex steroids and/or glucocorticoids). Using OCPI microscopy and a broad panel including both BAs and other polar steroids, we found a high degree of segregation between BA-sensitive and sulfated steroid–sensitive VSNs (Fig. 4). This indicates that BA responsiveness is not a general feature of all steroid-sensitive VSNs. It could be the case, however, that the BAs used in this panel are not the most potent members of the BA class. Future studies will be needed to determine the full breadth of VNO BA sensitivity, as there are dozens to hundreds of naturally occurring BA molecules that could have biological relevance in rodents (31).

In the initial experiments describing AOS BA sensitivity, a small number of AOB neurons (downstream of VSNs) demonstrated coactivation by BAs and sulfated glucocorticoids at 10 μM (4). AOB mitral cells innervate multiple AOB glomeruli and are capable of performing excitatory integration (22, 32), so it could have been the case that AOB mitral cells generated these response patterns through excitatory integration of glomerular input from BA- and glucocorticoid-specialist VSNs. Although we saw largely segregated VSN sensitivities to BAs and other steroids in this stimulus panel, we identified several VSN populations that were sensitive to both BAs and sulfated glucocorticoids (Fig. 4). This shows that some VSNs (and presumably the receptors they express) are capable of cross-class tuning upstream of the AOB circuit. We saw a very small number of VSNs that were coactivated by BAs and sulfated androgen, estrogen, and pregnanolone ligands, but given the relatively narrow focus of this panel, we cannot exclude that other non-BA steroid ligands may coactivate BA-sensitive VSNs. Overall, these studies show that the VSN sensitivity and selectivity for BAs is comparable to most of the other major known AOS ligand classes (6, 9, 16, 33, 34).

FLICCR-seq supports the identification of BA-sensitive VRs

An important step toward improving our understanding of mammalian BA chemosensation is to identify the peripheral receptors that endow animals with BA sensitivity. This continues to be a daunting task in the AOS, where limited tools have been available for receptor screening [although see (8, 11, 13)]. Rather than attempt to recreate heterologous expression systems (8, 35, 36), viral delivery methods (11), or transgenic mice for this purpose, we instead sought to develop a method that used primary VSNs in transgenic mice expressing GCaMP6f/s, drawing inspiration from (13). The major drawback in taking this approach is that, assuming proportional expression of the ~300 known VRs and FPRs (formyl peptide receptors) expressed by VSNs, less than 1% of the cells in the tissue express any given receptor (creating a “needle in a haystack” problem). We overcame this drawback by developing FLICCR-seq, which enabled us “find the needle in the haystack by making it glow.” In FLICCR-seq, we used clean micropipettes and an epifluorescence slice physiology microscope to pluck BA-responsive VSNs from VNO slices and process them for scRNAseq using established, commercially available protocols. FLICCR-seq samples were sequenced at relatively high coverage and with high quality (averaging over 90,000 mapped reads per sample), resulting in high-resolution coverage of expressed genes, including VRs, which are expressed at high relative levels in these cells (Fig. 5F). Analysis of FLICCR-seq VSN expression patterns produced strong matches with long-known signatures of V1R- and V2R-expressing cells (e.g., Gnai2 and Gnao1; Fig. 5), providing more evidence that FLICCR-seq produces interpretable data. In the context of the growing list of techniques for investigating VSN ligand-receptor interactions, the evidence indicates that FLICCR-seq has several benefits that make it an attractive approach in this difficult-to-study sensory system.

Despite its effectiveness in the context of this study, note a limitation of our initial implementation of FLICCR-seq, which relies on manual plucking of cells. Unlike droplet-based scRNAseq pipelines [reviewed in (37)], selecting cells manually from slices introduces the possibility that individual samples may contain the cell of interest plus some unintended material (e.g., debris from nearby dead cells). Although practice can reduce these occurrences, completely eliminating this material is nearly impossible. Rather than throw away samples in which a small amount of unintended material was visible, we instead focused on maintaining thorough annotation of each cell picking event, for example, by taking continuous screen captures of the procedure and making copious contemporaneous notes (fig. S6). Since the goal of FLICCR-seq was to identify the enrichment of receptors in BA-sensitive cells, so long as the identity (and dominant sensory receptor) of any unintentionally collected material was random, we would not expect that contaminating RNA sequencing counts would covary with annotated ligand sensitivities. Hence, in this context, there is little risk of false positives. Even still, improvements in the procedures for cell identification and clean sample picking will make FLICCR-seq even stronger and will support the discovery of more ligand-receptor matches in the future.

Identification of five V1Rs sensitive to four common BAs

The FLICCR-seq process allowed us to perform differential expression analysis on BA-sensitive and BA-insensitive VSNs (Figs. 5 and 6 and tables S1 and S2). Using multiple analytical strategies, five candidate BA-sensitive VRs emerged: Vmn1r16 (V1rc29), Vmn1r21 (V1rc28), Vmn1r70 (V1rl1), Vmn1r225 (V1re5), and Vmn1r228 (V1re3). Vmn1r16 and Vmn1r21 are both members of the V1rc clade on Mus musculus chromosome 6 (38) and were previously associated with in vivo activation by soiled bedding from multiple species (10). These results suggest that components of soiled bedding samples from rodent and nonrodent species activate BA-sensitive VRs. Differential expression analysis revealed enrichment of CA-sensitive cells in Vmn1r16-expressing VSN samples (Fig. 6 and tables S1 and S2) and DCA-sensitive cells in Vmn1r21-expressing samples. These two BA ligands are closely related (DCA is a gut microbe–dependent metabolite of CA) and are common to many species (18). Given that we observed CA- and DCA-selective VSNs in OCPI imaging experiments (Figs. 2 to 4), it seems likely that these two V1rc clade members are expressed by these sensitive, selective VSN populations. For these and other identified BA receptors, it is possible that the “best ligand” for the receptor may be outside the ligands tested here. However, given that in vivo BA exposure causes IEG induction in Vmn1r16 and Vmn1r21-expressing VSNs, it seems likely that these receptors play roles in sensing common environmental BAs. Similarly, differential FLICCR-seq analysis indicates that Vmn1r70 (V1rl1)–expressing VSNs are sensitive and selective for LCA, a potentially toxic natural metabolite of CDCA that is absent in rodents (Fig. 6 and tables S1 and S2) (18). OCPI experiments indicated the presence of VSNs that were sensitive and selective for 1 μM LCA, although this ligand may have caused some nonspecific cytotoxicity at 10 μM (Figs. 3 and 4 and figs. S3 and S4). To date, there is little known about V1rl1, located on M. musculus chromosome 7 adjacent to VRs assigned to the V1re clade (38). Our results suggest that V1rl1 may be a selective detector of LCA.

Differential gene expression analysis of Vmn1r225- and Vmn1r228-expressing samples indicated enrichment for multiple BA-sensitive VSNs (Fig. 6 and tables S1 and S2). Vmn1r225 (V1re5) and Vmn1r228 (V1re3) are both located on M. musculus chromosome 7 and have been associated with the V1re clade (38). These receptors were the most highly enriched in our FLICCR-seq dataset (Fig. 6), which was, by design, targeted at BA-sensitive VSNs. Given this experimental design, one would expect to encounter a higher number of VSNs expressing broadly tuned receptors than VSNs expressing narrowly tuned receptors. OCPI analysis, similarly, identified multiple populations of broadly tuned BA receptors (Fig. 4). These results suggest that Vmn1r225 and Vmn1r228 encode broad BA receptors. Previous studies indicated that Vmn1r225 and Vmn1r228 are activated by animal bedding (10). Earlier studies also indicated that closely related V1re clade members, Vmn1r226 (V1re2) and Vmn1r227 (V1re6), were sensitive to sulfated glucocorticoids (10). Given that OCPI experiments indicated that many broadly BA-tuned VSNs are also sensitive to 10 μM sulfated corticosterone (Fig. 4), it seems likely that these closely related receptors are capable of sensing steroids across the BA and glucocorticoid ligand classes.

Note that we encountered many samples in our FLICCR-seq dataset that expressed Vmn1r89 (V1rj2), a receptor that has been associated with sulfated estrogen and androgen sensing (Fig. 6) (10, 13). However, differential gene expression analysis did not indicate enrichment of BA-sensitive VSNs in Vmn1r89-expressing FLICCR-seq samples (tables S1 and S2). It is unclear why our FLICCR-seq dataset included so many Vmn1r89-expressing samples. A small number of sulfated estradiol- or testosterone-sensitive VSNs had some BA sensitivity (Fig. 4), but these data do not indicate a specific role for this receptor in BA detection.

Overall, the FLICCR-seq approach identified five strong BA receptor candidates, all of which were validated by duplex RNA in situ hybridization following in vivo BA exposure (Fig. 6). These receptors include members of V1R clades that have previously been associated with steroid sensation (e.g., the V1re clade) (10) and other clades that had not previously been associated with any specific ligand (e.g., V1rl and V1rc clades). The unbiased nature of the FLICCR-seq approach not only has the benefit of discovering receptor-ligand pairs that were not predicted by previous studies but also comes with limitations. For example, it is relatively laborious and expensive, which can impart restrictions on sample sizes. Similarly, using this approach, one cannot exclude the possibility that some BA-sensitive receptors were missed (e.g., lowly expressed mRNAs or receptors expressed in extremely rare populations). Even still, the benefits of this approach far outweigh the detractors and make this approach a compelling platform for future ligand-receptor discovery studies in the VNO and elsewhere in the body.

Together, the combination of OCPI and FLICCR-seq experiments produced a deep demonstration of BA sensitivity and selectivity and identified several V1R family members that serve as mouse BA receptors. The behavioral impacts of chemosensory BA detection remain unknown, but given the breadth of naturally occurring BAs and their associations with physiological states like species, sex, diet, and gut microbiome (18), it seems likely that BAs contribute to large repertoire of AOS-associated behaviors. The idea that BAs could play important roles in vertebrate social/reproductive behaviors is not new, as BAs are known pheromones in lampreys and other fishes (27). These studies indicate that these roles may also be played by BAs through the AOS in rodents. Overall, these results will support future studies into the molecular mechanisms of BA sensitivity, the network mechanisms that compare and combine BA information with other AOS ligands, and the impacts of BA chemosensation on mouse social behavior.



All animal procedures were performed in accordance with the Institutional Animal Care and Use Committee at the University of Texas Southwestern Medical Center and follow guidelines from the National Institutes of Health. Physiological experiments were performed with OMPtm4(cre)Mom/J knock-in mice (OMP-Cre mice; the Jackson Laboratory stock no. 006668) mated to either Gt(ROSA)26Sortm96(CAG-GCaMP6s)Hze/J mice (Ai96 mice; the Jackson Laboratory stock no. 024106) or CG-Igs7tm148.1(tet0-GCaMP6f,CAG-tTA2)Hze/J mice (Ai148 mice; the Jackson Laboratory stock no. 030328). These mice express the genetically encoded Ca2+ indicator GCaMP6s (OMP-Cre+/−, Ai96+/−) or GCaMP6f (OMP-Cre+/−, Ai148+/−) in VSNs, further referred to as OMPxAi96 mice or OMPxAi148 mice, respectively. Mice were between 6 and 15 weeks of age. The number, strain, and sex of the animals used in each experiment are described in corresponding figure legends. BaseScope experiments were performed with wild-type C57Bl/6J animals.

Solutions and stimulus presentation

All reagents were purchased from Sigma-Aldrich (St. Louis, MO, USA) unless otherwise specified. BA stimuli included CA, DCA, CDCA, and LCA. Sulfated steroids included epitestosterone-17 sulfate (A6940), 17α-estradiol-3-sulphate (E0893), 5α-pregnan-3β-ol-20-one sulfate (P3865), and corticosterone-21-sulfate (Q1570) purchased from Steraloids Inc. (Newport, RI, USA). Stock solutions (20 mM) of all BAs and sulfated steroids were prepared in methanol and diluted to their final concentration in Ringer’s solution containing 115 mM NaCl, 5 mM KCl, 2 mM CaCl2, 2 mM MgCl2, 25 mM NaHCO3, 10 mM Hepes, and 10 mM glucose. Artificial cerebrospinal fluid (aCSF) contained 125 mM NaCl, 2.5 mM KCl, 2 mM CaCl2, 1 mM MgCl2, 25 mM NaHCO3, 1.25 mM NaH2PO4, 25 mM glucose, 3 mM myoinositol, 2 mM sodium pyruvate, and 0.5 mM sodium ascorbate. All sulfated steroids were diluted to 10 μM (1:2000) for experiments. BAs were diluted to a range of concentrations from 0.1 to 10 μM. Control stimuli consisted of Ringer’s solution with 1:2000 methanol (the highest concentration of methanol in any individual stimulus). Stimuli were applied for 15 s using an air pressure–driven reservoir via a 16-in-1 multi-barrel “perfusion pencil” (AutoMate Scientific, Berkeley, CA, USA).

Volumetric VNO Ca2+ imaging

Following deep isofluorane anesthesia and rapid decapitation, VNOs were dissected, and the vomeronasal epithelium was carefully removed under a dissection microscope (Leica Microsystems, Buffalo Grove, IL, USA) similar to previous studies (19, 29). The vomeronasal epithelium was mounted onto nitrocellulose paper (Thermo Fisher Scientific, Atlanta, GA, USA) before placement into a custom imaging chamber. Volumetric Ca2+ imaging was performed using a custom OCPI microscope (30) with refinements described previously (29). GCaMP6s/f fluorescence was acquired using custom software that synchronized imaging with stimulus delivery via a randomized, interleaved stimulus delivery system (AutoMate Scientific). Images stacks containing 51 frames and spanning ~700 μm laterally, 250 to 400 μm axially, and ~150 μm in depth were taken once every 3 s (~0.33 Hz). Each individual stimulus was delivered for five consecutive stacks (~15 s) with at least 10 stacks between stimulus trials (≥30 s). Each analyzed experiment completed at least three complete randomized, interleaved stimulus blocks. Stimulation patterns are further described in Results.

Data analysis of volumetric VNO Ca2+ imaging

Data analysis was performed using custom MATLAB software similar to previous studies (19, 23). Briefly, image stacks underwent a rigid registration followed by a nonrigid warping. Then, ΔF/F, the relative change in GCaMP6s/f intensity, was calculated by subtracting the mean voxel intensity in three consecutive prestimulus stacks from the mean voxel intensity of three stacks during stimulus delivery and then dividing by the value of the mean prestimulus intensity. Volumetric ROIs were manually drawn around the cell bodies of well-registered VSNs that reliably responded to stimulation as well as spontaneously active neurons. Following the volumetric ROI selection, the mean voxel intensity for each ROI was calculated for every image stack in the experiment (~1200 stacks), generating a matrix of fluorescence intensity. The across-trial mean ΔF/F for all ROIs and stimulus applications was calculated using this matrix. Stimulus responsiveness was assessed by comparing the across-trial stimulus-evoked ΔF/F of each ROI to the across-trial Ringer’s control stimulus (to account for potential effects of valve switching). The nonparametric Wilcoxon-Mann-Whitney test was used to assign a P value associated with the stimulus response, and ROIs demonstrating a positive ΔF/F response and a P < 0.05 to a stimulus were considered responsive. For clustering analysis, the normalized stimulus-evoked ΔF/F of all ROIs were clustered using a bootstrapping method based on mean shift (22, 39). Cluster separation was assessed using linear discriminant analysis (LDA). Within each experiment, each group of clustered responses was projected along the calculated pairwise LDA eigenvector. Discriminability index (d′) was calculated on the basis of the LDA projections of each group, resulting in a matrix of pairwise cluster separation (Figs. 2E, 3E, and 4E). A discriminability index score exceeding 1.64 (when the z-score equivalent P value is equal to 0.05) was used to confirm satisfactory cluster separation.

Ex vivo Ca2+ imaging and data analysis

Glomerular imaging was performed using a VNO-AOB ex vivo preparation (23, 40). Briefly, after deep isofluorane anesthesia, mice were decapitated and placed into ice-cold aCSF. The anterior skull including the snout was dissected from the skull and then hemisphered along the midline. The single hemisphere of the mouse snout including the olfactory bulb was adhered to a small plastic plank and submerged in rapidly circulating aCSF in a custom dissection chamber. Under a dissection microscope, the septal cartilage and bone overlying the vomeronasal nerve was carefully removed, the AOB separated from the frontal neocortex, and a 0.0045″ OD (outer diameter) polyimide cannula (A-M Systems, Carlsborg, WA, USA) was then placed into the VNO for ligand delivery. Volumetric images of the AOB were then acquired using OCPI microscopy using the same stimulus protocols described for VNO imaging. Maps of volumetric glomerular activation were generated as described in (23) (fig. S5).

VNO slicing and single VSN collection

For FLICCR-seq experiments, acute VNO slices from OMPxAi96 mice were prepared on a vibrating microtome (Model VT1200S, Leica) and placed into a modified slice chamber (Warner Instruments). An upright electrophysiological rig based on the Nikon FN1 housing was modified to accommodate a custom computer-controlled stimuli delivery system. At least three randomized, interleaved trials of BA stimuli (Ringer’s, 1 μM CA, 1 μM DCA, 1 μM CDCA, and 1 μM LCA) were applied to the tissues. Because of technical limitations, each VNO slice was exposed to two of the four tested BA ligands (either 1 μM CA and 1 μM DCA or 1 μM CDCA and 1 μM LCA; Figs. 5F and 6A). A heat map representation of the change in fluorescence intensity (ΔF/F) was generated in near real time by a customized MATLAB script to identify reliably responsive VSNs in the tissue. Responsive VSNs showed stimulus-evoked GCaMP6s increases to at least three randomly presented trials of the test BA and did not respond to Ringer’s control stimulation.

A modified patch pipette (tip diameter, 1.5 to 2 μm) was placed near responsive neurons by comparing the ΔF/F image to infrared differential interference contrast images. Verification of the pipette placement near the intended VSN was performed by pulsing the tissue with BA stimuli, while the pipette was within 5 μm of the cell of interest. Once the location of the BA-responsive VSN was confirmed, the VSN soma was gently aspirated into the pipette tip and then placed into RNA preservation buffer for later scRNAseq. Throughout the single cell selection procedure, images and movies were taken of the cellular and pipette positions, and careful notes were taken about the responsiveness of the VSN, whether any unintended material (e.g., an adjacent dead cell on the outside of the pipette) was visible. Successful capture of the intended VSN was assessed by visually observing the fluorescent, stimulus-responsive cell within the patch pipette before retrieval and then observing the absence of the intended stimulus-responsive cell in the tissue after pipette retrieval.

mRNA from each individual VSN sample was converted to complementary DNA (cDNA) using the SMART-Seq v4 Ultra Low Input RNA Kit (Takara Bio USA Inc.). The quality of each individual cDNA sample was assessed using an Agilent TapeStation, and samples lacking cDNA fragments spanning the 500 to 2500 base pair (bp) size were excluded from further processing and analysis. cDNA libraries were then prepared using the Nextera XT DNA Library Prep Kit (Illumina, CA, USA) and evaluated again for fragment size using the Agilent TapeStation. Each individual sample cDNA library was barcoded using a Nextera XT Index kit (Illumina) and diluted to approximately match concentrations across samples. Indexed sample libraries were pooled and sequenced on the MiSeq platform (150 bp, paired end, ~12 million reads per pooled library) by the staff of the McDermott Next Generation Sequencing Core at University of Texas Southwestern Medical Center.

scRNAseq analysis

Individual fastq files for each VSN sample were initially generated from raw data by bioinformaticians in the McDermott Next Generation Sequencing Core. Samples were trimmed via Trimmomatic (41) and aligned to the mm10 M. musculus reference genome using the STAR aligner (42). A total of 158 samples were analyzed, including 145 single VSN samples, seven bulk VNO tissue samples, and six kit controls. The average number of mapped reads per sample was 90,636 (median, 83,955). A matrix of raw counts per gene (using HTseq, ENSEMBL version 90) for each sample was generated and imported into Seurat v 2.3.4 (24) for subsequent analysis. Briefly, the raw reads were normalized using the “lognormalize” subroutine, the top differentially expressed genes were selected using the “FindVariableGenes” function, and the resulting normalized data scaled using the “ScaleData” function. Principal components analysis (PCA) was run on the variable genes, and approximately 10 PCs were selected that contained the most useful gene expression eigenvectors. Clustering was performed using the default routines of the “FindClusters” function (“resolution” value 1.0), and data were prepared for visualization using UMAP (uniform manifold approximation and projection) (43). Differential gene expression analysis comparing BA-responsive and BA-unresponsive populations was conducted using the Wilcoxon-Mann-Whitney test (via the “FindMarkers” function in Seurat). For this comparison, the log fold change threshold was set at 0.5, and the minimum percentage of expressing samples set to 0.01 (i.e., at least one sample). The differentially expressed genes with unadjusted differential expression P < 0.001 were sorted by log fold change (tables S1 and S2).

For additional analysis of gene expression patterns within the vomeronasal gene family, a subset of the log-normalized data including gene expression information for VRs, FPRs, and selected housekeeping and VNO-associated genes (including Actb, Gapdh, Gnai2, Gnao1, Gnaq, Gnas, Obp2a, Obp2b, Omp, Osbp, Osbp2, Slc9a3r1, and Trpc2) were subjected to an independent cluster analysis using previously reported routines based on the “mean-shift” method (22, 39). Clusters identified using this alternative method were organized for display purposes, and the binomial distribution were used to quantify the likelihood of specific clusters containing the observed number of BA-responsive VSNs.

In vivo BA exposure and in situ hybridization

In situ hybridization was performed on VNO slices obtained from 24 C57Bl/6J mice (12 males and 12 females). Briefly, mice were lightly anesthetized via isofluorane, and 20 μl of one of the follow stimuli was pipetted directly onto the mouse snout (~10 μl per nostril); Ringer’s, 1 mM BA mix, 100 μM BA mix, or female mouse urine (not shown). Mice were returned to clean cages and perfused with 4% paraformaldehyde 30 min following exposure. The VNOs were extracted and later sectioned using a cryostat (Model CM-1950, Leica, Buffalo Grove, IL) and mounted onto slides. Custom BaseScope probes for Vmn1r16, Vmn1r21, Vmn1r70, Vmn1r225, Vmn1r228, and Vmn2r54 were generated by a commercial vendor (Advanced Cell Diagnostics, Newark, CA). Because of the high degree of homology between receptors of the same clade, the two V1rc clade probes may have cross-detected other V1rcs. Specifically, the probe for Vmn1r16 could cross-detect Vmn1r22, Vmn1r23, Vmn1r24, and Vmn1r29, and the probe for Vmn1r21 could cross-detect Vmn1r27 and Vmn1r28. BaseScope probes for the IEGs Egr1 and Fos were existing commercial products that we purchased without modification. The VR probes were labeled with one colorimetric label, and the two IEG probes were both marked with a second colorimetric label using a BaseScope Duplex Detection Reagent Kit (Advanced Cell Diagnostics, Newark, CA) following the manufacturer’s instructions. The samples were counterstained with dilute hematoxylin and eosin and slides imaged using a Nanozoomer HT 2.0 slide scanner (Hamamatsu, Japan) maintained by the University of Texas Southwestern Whole Brain Microscopy Facility.

In situ hybridization quantification

Nanozoomer images were exported to high-resolution color tiff images and then processed using publicly available image classification software, ilastik. Briefly, the “Pixel Classification” plugin for ilastik, which implements voxel-wise, supervised random forest machine learning, was used to produce probability map images of the positions of VR- and Egr1-positive puncta. Representative blocks of tissue sections for each VR probe were used as model training datasets, and the same model was applied to 1067 VNO tissue sections across 116 tissue blocks derived from 24 experimental animals. The probability map images were then used as the input images for CellProfiler, which was used to segment IEG- and VR-positive areas. Last, the output images from CellProfiler were assembled via custom a MATLAB user interface that allowed users to draw and annotate ROIs for each imaged vomeronasal epithelium (e.g., sex and treatment condition). MATLAB scripts were used to assemble the data and perform statistical comparisons. Colocalization of IEG and VR probes was assessed by identifying the presence of one or more IEG puncta within 6 μm of the center of a VR-positive area (i.e., within the soma radius of a VSN expressing that VR). The block-wise design of our section mounting strategy allowed unambiguous assignment of tissues to their proper treatment condition and sex but did not allow us to unambiguously assign each section to a specific animal. We pooled sections by condition and sex and then performed statistical comparisons using the Kruskal-Wallis test, followed by post hoc multiple comparisons tests as appropriate.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank N. Browder and C. Nielson for technical support and G. Konopka and F. Ayhan for critical feedback on the manuscript. Funding: This work was supported by the National Institute on Deafness and Other Communication Disorders and National Institute of Neurological Disorders and Stroke of the U.S. NIH [grants R01DC015784 (to J.P.M.), R01DC017985 (to J.P.M.), R21NS104826 (to J.P.M.), and F31DC017661 (to W.M.W.)] and the Welch Foundation [grant number I-1934-20170325 (to J.P.M.)]. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. Author contributions: W.M.W., J.C., X.Z., W.I.D., L.L.M., and L.G. designed and performed experiments. W.M.W., J.C., and J.P.M. developed FLICCR-seq. W.M.W., J.C., W.I.D., and J.P.M. analyzed data. W.M.W., J.C., and J.P.M. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article