Abstract
Summary
We have generated single cell transcriptomic atlases of vomeronasal organs (VNO) from juvenile and adult mice. Combined with spatial molecular imaging, we uncover a distinct, previously unidentified class of cells that express the vomeronasal receptors and a population of canonical olfactory sensory neurons in the VNO. High resolution trajectory and cluster analyses reveal the lineage relationship, spatial distribution of cell types, and a putative cascade of molecular events that specify the V1r, V2r, and OR lineages from a common stem cell population. The expression of vomeronasal and olfactory receptors follow power law distributions, but there is high variability in average expression levels between individual receptor and cell types. Substantial co-expression is found between receptors across clades, from different classes, and between olfactory and vomeronasal receptors, with nearly half from pairs located on the same chromosome. Interestingly, the expression of V2r, but not V1r, genes is associated with various transcription factors, suggesting distinct mechanisms of receptor choice associated with the two cell types. We identify association between transcription factors, surface axon guidance molecules, and individual VRs, thereby uncovering a molecular code that guides the specification of the vomeronasal circuitry. Our study provides a wealth of data on the development and organization of the accessory olfactory system at both cellular and molecular levels to enable a deeper understanding of vomeronasal system function.
Introduction
In many terrestrial species, the vomeronasal organ (VNO) is dedicated to the detection of inter- and intra-species chemosensory cues1–5. Detection of these cues triggers innate neuroendocrine responses and elicits stereotypic social and reproductive behaviors1,2,6–16. The VNO shares a developmental origin with the main olfactory epithelium (MOE), which detects the odor world at large and allows associative learning to take place. Both develop from the olfactory placode during early embryogenesis17,18, but the two sensory organs follow different developmental trajectories to establish distinct characteristics in morphology, cellular composition, and molecular features. Single cell atlases of the MOE have been generated to reveal astonishing details in its molecular composition and developmental trajectories19–23. Transcriptomic data of the VNO is not extensive24,25, but single cell analysis have already provided critical information about VNO development 26,27. In this study, we generate single cell atlases of developing and adult mouse VNO neuroepithelia to answer fundamental questions about the VNO.
The vomeronasal sensory neurons (VSNs) express three families of G-protein coupled receptors, the V1rs, V2rs, and the formyl peptide receptors (Fprs)28–34. With more than 400 members, the vomeronasal receptors are among the fastest evolving genes35–40. Signaling pathways in the VNO include the Gi2 and Go proteins, and a combination of Trpc2, Girk1, Sk3, and Tmem16a ion channels, to transduce activation of the vomeronasal receptors41–61. The spatially segregated expression of Gi2 and Go, as well as the VR genes suggests two major classes of neurons. Our study reveals new classes of sensory neurons, including a class of canonical olfactory sensory neurons (OSNs) in the VNO. We further determine the developmental trajectories of the separate lineages and the transcriptional events that specify these lineages.
Pheromones are highly specific in activating the VSNs5,62–67. Previous studies have shown that VSNs residing in the apical layer of the VNO express V1R (Vmn1r) genes in a monoallelic manner65,68,69, whereas the basal VSNs express one or two of the broadly expressed C clade of V2R (Vmn2r) genes, plus a unique V2R gene belonging to another clade28–33,70–72. Although these results suggest that receptor expression in the apical VNO conforms to the “one neuron one receptor” pattern as found in the MOE, the mechanisms that control receptor expression are unknown. Here, we find substantial co-expression of vomeronasal receptors, and of vomeronasal and odorant receptors. Moreover, our analyses indicate that selection of V1R expression likely results from stochastic regulation as in the OSNs, but V2R expression likely result from deterministic regulation.
Finally, we address the molecular underpinning of how VSNs establish anatomical connections to transmit sensory information. VSNs expressing the same receptor project to dozens of glomeruli in the AOB69,73. Individual stimuli activate broad areas in the AOB74. Moreover, the dendrites of the mitral cells in the AOB innervate multiple glomeruli75–77. This multi-glomerular innervation pattern is in stark contrast with the main olfactory system, where olfactory sensory neurons (OSNs) expressing the same odorant receptor converge their axons into mostly a single glomerulus in each hemisphere of the main olfactory bulb (MOB)78. The anatomical arrangement in the AOB has strong implications as to how species-specific cues are encoded and how the information is processed. In the MOB, when the convergent glomerular innervation is experimentally perturbed to become divergent, it does not affect detection or discrimination of odorants but diminishes behavioral responses to innately recognized odors79,80. Thus, stereotypic projection patterns provide a basis for genetically specified connections in the neural circuitry to enable innate behaviors. Consistent with this notion, it has been shown that mitral cell dendrites innervate glomeruli containing the same vomeronasal receptor (VR) type such that the divergent projection pattern of the VSNs is rendered convergent by the mitral cells75. This homotypic convergence suggests that rather than using the spatial position of the glomeruli, the connection between VSNs and mitral cells in the AOB may rely on molecular cues to enable innate, stereotypical responses across different individuals. To a lesser extent, heterotypic convergence, i.e., axons expressing different receptors innervating the same glomeruli, is also observed76. In either case, expression of the molecular cues is likely genetically specified and tied to individual VRs, but little is known about how this specification is determined. Our analyses revealed the stereotypic association between transcription factors, axon guidance molecules, and the VRs to suggest a molecular code for circuit specification. Whereas this manuscript highlights some of the main discoveries, much detailed analyses can be found in the dataset hosted online for readers to browse.
Results
Cell types in the VNO
We dissected mouse VNOs from postnatal day 14 (P14) juveniles and P56 adults. Cells were dissociated in the presence of actinomycin D to prevent procedure-induced transcription. From four adult (P56) and four juvenile (P14) mice (equal representation of sexes) we obtained sequence reads from 34,519 cells. The samples and replicates were integrated for cell clustering. In two-dimensional UMAP space, 18 cell clusters can be clearly identified (Fig. 1A). These clusters were curated using known cell markers (Fig. 1B). There was no obvious difference in the presence of cell clusters between juvenile and adult VNOs (Fig. 1C) or between male and female sexes (Fig. 1D), though there were differences in gene expression profiles between the ages (Fig. S1) and sexes (Fig. S2), the list of significantly differentially expressed genes did not appear to be influential for the neuronal lineage and cell type specification, or related to cell adhesion molecules, which were the main focuses of this study.
Clustering confirmed all previously identified cell types but also revealed some surprises. The largest portion of cells belonged to the neuronal lineage, including the globose basal cells (GBC), immediate neuronal progenitors (INPs), the immature and mature VSNs, and a population of cells that do not co-cluster with either VSN type (see below). There were substantial numbers of sustentacular cells (SCs), horizontal basal cells (HBCs), and ms4-expressing microvillus cells (MVs). Cells engaged in adaptive immune responses, including microglia and T-cells, were also detected. A population of Fpr-1 expressing cells that were distinctive from the VSNs expressing the Fpr family of genes formed a separate class. These were likely resident cells mediating innate immune responses. We also identified the olfactory ensheathing cells (OECs) and a population of lamina propria (LP) cells, which share molecular characteristics with what we have found in the MOE81.
To obtain the spatial location of the various cell types, we selected 100 target genes based on the scRNA-seq results. Using probes for these genes, we used the Molecular Cartography® platform to perform spatial molecular imaging (Fig. 1E). Based on molecular clouds and DAPI nuclei staining, we segmented the cells and quantified gene expression profiles to cluster the cells. We then map individual cell clusters onto their spatial locations in VNO slices. Unlike previous studies that relied on a few markers to identify cell types, our approach relied on the spatial transcriptome to calculate the probability that a cell belongs to a specific class. This analysis revealed that the VSNs and supporting cells are located in the pseudostratified neuroepithelium (Fig. 1F). The LP cells are located along the lamina propria underlying the neuroepithelium as found in the MOE (Figs. 1F and 1G). Surprisingly, however, there are few HBCs located along the basal lamina, in direct contrast to the MOE (Figs. 1F and ="fig"> 1G). Most of the HBCs are found in the non-neuronal epithelium surrounding the blood vessel, with few near the marginal zone. The marginal zone is thought to be the neurogenic regionError! Hyperlink reference not valid.. We quantified the distribution of in various regions of the VNO neuroepithelia (Fig. S3A) and found significantly more GBCs, INPs, and immature VSNs in marginal zone than in the main zone (Figs. S3B, and S3C). A previous study suggested neurogenic activity in the medial zone82, but we did not find evidence in support of this conclusion, which was largely based on BrdU staining of mitotic cells without lineage-specific information. Based on our transcriptomic analysis, we conclude that neurogenic activity is restricted to the marginal zone.
Novel classes of sensory neurons in the VNO
To better understand the developmental trajectory of the VSNs, we segregated the cells in the neuronal lineage from the whole dataset (Fig. 2A). The neuronal lineage consists of the GBCs, INPs, immature neurons as determined by the expression of Gap43 and Stmn2 genes (Figs. S4A and S4B), and the mature VSNs. Cells expressing Xist, which is expressed in female cells, were intermingled with those from males (Fig. S4C). This observation is consistent with our previous studies using bulk sequencing indicating that the cell types are not sexually dimorphic25. It is also consistent with physiological responses of the VSNs to various stimuli5,64,66,83,84. For further analyses, therefore, we did not segregate the samples according to sex.
The VSNs are clearly separated into the V1R and V2R clusters as distinguished by the expression of Gnai2 and Gnao1, respectively (Fig. 2B). Surprisingly, we observed that a portion of the V2R cells also expressed Gnai2. These cells were primarily from adult, but not juvenile male mice (Fig. S4D). Past studies based on Gnai2 and Gnao1 expression would have identified this group of cells as V1R VSNs, but the transcriptome-based classification places them as V2R VSNs. The signaling mechanism of this group of cells may be different from the canonical V2R VSNs.
We found a major group of cells that expressed the prototypical VSN markers but formed a cluster distinct from the mature V1R and V2R cells (Figs. 1A and 2A). Within the cluster, there was an apparent segregation among the cells into V1R and V2R groups based on marker gene expression (Fig. 2B). A small group of OR-expressing neurons also belonged to this cluster. These cells expressed fewer overall genes and with lower total counts when compared with the mature V1R and V2R cells (Figs. S4E and S4F). The lower ribosomal gene expression suggests that these neurons are less active in protein translation (Fig. S4G). This group of cells are scattered within the epithelium (Fig. 2E).
There are 503 differentially expressed genes in this cluster when compared with other mature neurons (padj < 0.001; FC > 1.5; Figs. 2F, S4H). Gene ontology (GO) analysis reveals that differentially expressed genes enriched in this group are associated with odorant binding proteins (Fig. 2G). Correspondingly, mucin 2 (Muc2), odorant binding protein 2a (Obp2a), Obp2b, and lipocalin 3 (Lcn3) genes, usually enriched in non-neuronal supporting or secretory gland cells, are expressed at higher levels in these neurons than the canonical VSNs (Figs. 2H, S5A-D). They also have higher expression of Wnk1, which is involved in ion transport, the lncRNA Neat1, the purinergic receptor P2ry14, the zinc finger protein Zfp738, and the centrosome and spindle pole associated Cspp1 (Fig. S5E-I). On the other hand, these cells exhibit lower expression of Omp and JunD (Figs. 2F, S5M-N). The lower level of Omp suggests that these cells do not have the full characteristics of mature VSNs (mVSNs), but they are also distinct from the immature VSNs (iVSNs) as they do not express higher levels of immature markers such as Gap43 and Stmn2 (Figs. S4A and S4B). The lower expression of Ndua2, Ndua7, Cox6a1, and Cox7a2, which are involved in mitochondrial activities, suggest that these cells are less metabolically active than the canonical VSNs (Figs. 2F, S5O-R). To rule out the possibility that these cells were immature or senescent, we performed a pseudotime analysis on the neuronal lineage and found the cluster to have similar pseudotime values as the mature V1R and V2R linages (Fig. S6A and 6B). Taken together, the gene expression profile suggests that this new class of neurons may not only sense environmental stimulus, but also may provide proteins to facilitate the clearing of the chemicals upon stimulation. We, therefore, name these cells as putative secretory VSNs (sVSNs).
We also identify a distinct set of cells expressing the odorant receptors (ORs) that comprise ∼2.3% of the total neurons (Fig. 2A). A previous study has found the expression of ORs in neurons that display some typical molecular features of VSNs, while another study detecting ORs in the VNO suggested the presence of non-canonical OSNs 85,86 The OR expressing cells were shown to project to the AOB, but it was not clear how prevalent OR expression was in the VNO, nor whether the cells were VSNs or OSNs. While some OR expressing cells cluster with the mature V1R or V2R neuronal lineages, a majority of these cells lacks V1R or V2R markers and forms a cluster distinct from the V1R and V2R VSNs (Fig. 2C). These cells express Gnal and Cnga2, which are the canonical markers of OSNs in the main olfactory epithelium. Differential gene expression analysis (Fig. S6C) reveals an enrichment for multiple GO terms related to cilium (Fig. S6D), consistent with the ciliated nature of OSNs. We thus mark the cells as canonical OSNs. Spatial mapping reveals that the OSNs are mainly in the neuroepithelium, with some cells concentrated in the marginal zone (Fig. 2D).
Developmental trajectories of the neuronal lineage
The vomeronasal organ develops from the olfactory pit during the embryonic period and continues to develop into postnatal stages87–89. Neurons regenerate in adult animals90,91. Cell types in the vomeronasal lineage have been shown to be specified by BMP and Notch signaling26,92, and coordinated by transcription factors (TFs) including Bcl11b/Ctip2, C/EBPγ, ATF5, Gli3, Meis2, and Tfap2e93–98. Recent scRNA-seq studies of the VNO have helped identify Notch signaling as a specifier of the apical and basal lineages and have provided insight into the distinct transcriptional profiles of the basal and apical program26,27. Despite these advances in understanding the role of individual genes in VNO development, the transcriptional program that specifies the lineages is not known.
To explore VSN development, we performed pseudotime inference analysis of the V1R and V2R lineages for P14 and P56 mice using Slingshot (Fig. 3A)99. We set GBCs as the starting cluster and mVSNs as terminal clusters. A minimum spanning tree through the centroids of each cluster was calculated using the first fifty principal components and fit with a smooth representation to assign pseudotime values along the principal curve of each lineage. Cell density plots for both the V1R and V2R lineage reveal a higher portion of immature VSNs at P14 than at P56. For the mature VSNs, the P14 samples have peaks at an earlier pseudotime than the P56 mice (Fig. 3B). This result indicates the VSNs in juvenile mice are developmentally less mature than their counterparts in adults, but these differences do not distinguish them in obvious ways (Fig. 1C).
On the other hand, we have identified dynamic changes in transcriptomes associated with the V1R and V2R lineages. There were 2,037 significantly differentially expressed genes between V1R and V2R lineages (Fig. 3C). While V1R and V2R lineages shared some of the genes in the early stages of development, most show distinct expression patterns between the two. To determine the transcriptional program associated with the lineages, we examined the expression sequence of individual transcription factors and signaling molecules from the gene list.
The transcription factor Ascl1 is expressed by a subset of the GBCs that appears to be the earliest in developmental stage whereas Sox2 is expressed by a broader set of GBCs (Fig 3D-F). The early INPs are defined by the expression of Sp8, Neurog1 and NeuroD1. Neurog1 is mostly restricted in the early INPs whereas NeuroD1 and Sp8 are also expressed by the late INPs (Fig. 3G-I, S7A-B). NeuroD1 is expressed by iVSNs and iOSNs, but Sp8 is only expressed by the iVSNs.
To obtain a refined view of the transition from INPs to the immature sensory neurons, we took the subset and re-clustered them (Fig. 3J). The late INPs are separated into two clusters that share marker genes with the V1R and V2R iVSNs, respectively, indicating that commitment to the two lineages begins at the late INP stage. Consistent with previous findings, the homeobox protein Meis2 was specific to the V1R lineage (Fig. 3K). Concomitant with Meis2, there are a number of other genes expressed by the late INPs committed to the V1R fate, including secretin (Sct), Foxj1 target gene Fam183b, the microRNA Mir100hg, acyl-CoA dehydrogenase Acad10, and Keratin7 (Krt7; Fig. 3L-O, S7C). Notably, Bcl11b is exclusively absent from INPs committed to the V1R fate. This observation is consistent with the observation by Katreddi and colleagues that Bcl11b is lost in basal INPs in Notch knockout26. Together with the observation that loss of Bcl11b results in increased number of V1R VSNs93, these results indicate that transient downregulation of Bcl11b is required for the V1R lineage (Fig. 3U).
The transcription factor Tfap2e, which is required to maintain the V2R VSNs94, is expressed by the iVSNs but not by the late INPs committing to the V2R fate (Fig. 3P). Sp8 is the only transcription factor specifically expressed in the late INPs committed to the V2R but not the V1R fate (Fig. 3I). Emx2 is expressed by all neuronal lineage cells (Fig. S7E). Krt8 is also found throughout the V2R lineage, but its expression is diminished in the V1R cells (Fig. S7D).
The OSN lineage is marked by the expression of Olig2, Fezf1, Tshz2, and Nfib (Fig. 3Q-T). The expression of Olig2 and Fezf1 is exclusive to the OSN fate. Tshz2 is expressed at a late stage of the iVSNs, but not in the INPs. There are multiple genes that may not directly be engaged in cell fate determination but are clearly markers of the cell types (Fig. S7F-R). Although Notch1 and Dll4 are not identified as significantly differentially expressed, feature plots show clear distinction in their expression in the V2R and V1R lineage, respectively, as an earlier study showed26.
Based on these patterns, we propose a model of molecular cascades that specify the neuronal lineages in the VNO (Fig. 3V). Ascl1 and Sox2 specify the GBCs. The down regulation of Ascl1, Sox2, and subsequent upregulation of Neurog1, NeuroD1, and Sox8 commits the cells to the early INPs. The VSN and OSN lineages diverge after the early INP stage. The downregulation of Sp8, Nfib, and Bcl11b and the expression of Meis2 promote the V1R fate. The downregulation of Sp8 and the expression of Fezf1, Tshz2, and Olig2 are required for commitment to the OSN lineage. Downregulation of Sp8 is not required for the V2R lineage, which begin expressing Tfap2e. NeuroD1 is expressed in all INPs. These patterns of expression suggest that both the OSN and V1R lineages required the expression of specific transcription factors. The V2R lineage, on the other hand, appears to rely on factors inherited from the early INPs. This suggests the possibility that the V2R lineage is a default path for the VSNs.
Receptor expression in the VSNs
We quantified the expression of V1Rs, V2Rs, and ORs to gain insights into how chemosensory cues may be encoded by the VNO. For comparison, we also included OR expression from the MOE81. Within each class of receptors, the probability of expression of a gene follows a power law distribution except for the lower ranked genes (Fig. S8A). The sharp deviation from the power law curve for the lower ranked receptors is likely from technical dropout of genes expressed at low levels as they cannot be effectively captured by the scRNA-seq platforms. We plotted the relationship between total reads and the number of cells expressing a given receptor, and the average reads per cell for the receptors (Fig. 4A-H). We observed weak correlations between the total reads and the cell number expressing a given V1R or a V2R (Fig. 4A and B). This non-uniform distribution of VR genes is consistent with our observation from bulk sequencing results25.
Most V1Rs and V2Rs are expressed by less than 1500 cells (out of 34,519). Most VRs are expressed at less than 100 counts per cell (Fig. 4E and F). Several V1Rs, including Vmn1r184, Vmn1r89, Vmn1r196, Vmn1r43, and Vmn1r37, are highly expressed in individual cells (Fig. 4A). Vmn1r196, Vmn1r43, and Vmn1r37 are also expressed at the highest level per cell. Some others, including Vmn1r183, Vmn1r81, and Vmn1r13, are expressed in large numbers of cells but have low expression in individual cells. Notably, the two highest expressed receptors recognize female pheromone cues. Vmn1r89, also known as V1rj2, is one of the receptors that detects female estrus signals. Vmn1r184 is one of the receptors for the female identity pheromone4,63,100. The functions of Vmn1r196, Vmn1r43, and Vmn1r37, however, are unknown. We did not detect the expression of 45 V1Rs, which may be the result of technical dropout (Fig. 4E).
We do not observe obvious correlations between expression level and chromosomal locations. Both Vmn1r89 and Vmn1r184 are located on Chromosome 7, in a region enriched in V1R genes. Similarly, Vmn1r37 and Vmn1r43 are in a V1R-rich region on Chr. 6. However, Vmn1r183, Vmn1r13, and Vmn1r81, which are in the same clusters, are expressed by many cells but at some of the lowest levels.
All V2R genes are detected in the VNO (Fig. 4B and F). Vmn2r53, which has been shown to mediate intermale aggression through a dedicated circuit, has the highest level of expression and is expressed by the second most cells101. Among the highly expressed V2Rs, Vmn2r1 and Vmn2r7 are co-receptors for other V2Rs. Vmn2r59 has been shown to detect predator signals4. Vmn2r115, a receptor for ESP22 that is secreted by juveniles102, is expressed by the highest number of cells. However, Vmn2r116 (V2rp5), which recognizes ESP1 and induces lordosis behavior in females, is a close homolog of Vmn2r115 but not among the highly expressed genes102,103. Notably, Vmn2r114, close homolog of Vmn2r115 and Vmn2r116, is also expressed by large numbers of cells. Vmn2r88, a hemoglobin receptor104, was not identified as a highly expressed gene. The functions of other highly expressed receptors are not known.
In contrast to the VR genes, total counts for ORs in the VNO exhibit a tight relationship with the number of cells (Fig. 4C and G). Except for Olfr124, most of the OR genes align well with the linear regression curve. This relationship is different from the VR genes and is also distinct from the single cell data from the MOE, which exhibits a similar distribution as the VRs81 (Fig. 4D and H). Out of the 686 OR genes detected in the VNO, only 80 are expressed by more than 10 copies per cells, indicating that a majority of the ORs do not contribute to meaningful signaling of chemosensory cues.
To comprehensively survey receptor expression, we also included VR and OR pseudogenes in our analysis (Fig. S8B-G). We did not detect a significant expression of pseudogene V1Rs (Figs. S8B and S8E), but Vmn2r-ps87 has the highest count/cell in the V2R population (Fig. S8F). Olfr709-ps1, and Olfr1372-ps1 are the two highest expressed genes in terms of total count and total number of cells in the VNO (Fig. S8D).
We next examined transcription factors associated with individual receptor types. In the MOE, ORs are monoallelically expressed105. The unique expression of an OR gene is mediated by chromosomal repression and de-repression, coordinated by transcription factors and genes involved in epigenetic modification106–111. Monoallelic V1R gene expression is also observed69. Epigenetic modification takes place at V1R gene clusters112,113, but the repression of receptor genes appears to permit transcriptional stability rather than receptor choice114. How VR genes are selected is not known. To explore our dataset for clues of transcriptional activities associated with VR expression, we plotted the cross-correlation between VRs and their TF profiles. We observed correlations among receptors (Figs. 4I and 4J). We did not find an obvious association between TF profiles and VR sequence similarity of pairs (Fig. S6A).
There is a high correlation among a large portion of V1Rs (Fig. 4I). By analyzing the correlation between individual TFs with the VRs, we found that V1R expression is overwhelmingly associated with Meis2 (Fig. S9B). The few receptor types that do not show high correlation with Meis2 are associated with Egr1 or c-Fos. It is not clear whether these prototypical immediate early genes are involved in specifying receptor choice. This result indicates that once the V1R lineage is specified by Meis2, receptor choice is not determined by specific combinations of transcription factors. This scenario is similar to receptor choice by the OSNs in the MOE.
Different from the V1Rs, we observed islands of high similarity of TF expression among the V2Rs, indicating that receptors in these islands share the same set of TFs (Fig. 4J). We identify correlation between individual TFs with V2Rs (Fig. S9C). Whereas Tfap2e is involved in specifying the V2R fate, it is only associated with the expression of a subset of V2Rs. Pou2f1 and Atf5 are more strongly associated with other subsets of receptors. Unlike the V1Rs, there is a disparate set of TFs associated with the V2Rs, suggesting that the V2R choice may be determined by combinations of TFs. Based on the prevalence of individual transcription factors in association with the VRs, we propose a model of transcriptional cascade that may be involved in receptor choice (Fig. 4K). For the V2Rs, Tfap2e+ cells can be further specified by Ikzf4, Tcf4, and Trps1. In Ikzf4 negative cells, Batf3, Atf5, Pbx2, and Pou2f1 can specify receptor types, respectively. In the Tfap2e-cells, receptors can be specified by Pou2f1, Rlf, and Batf3.
Co-expression of chemosensory receptors
We next investigated co-expression of receptors in individual cells. Since we applied SoupX to limit ambient RNA contamination and Scrublet to remove doublet cells, we set a stringent criterion in counting receptors expressed by single cells115. V1R and V2R genes on average constitute ∼2% of total reads per cell, and the ORs constitute less than 1% of the total reads per cell (Fig. S10A-C). On average, the V2Rs have significantly more than one receptor per cell (Fig. S10B). Using Shannon Index to measure uniqueness of receptor expression in individual cells, we found that the mature VSN and OSNs have relatively high specificity, but many cells show significantly higher index values, indicating that they expressed multiple receptors (Fig. 5A). In support of this observation, there are significant representations of the second and third highest expressed receptors in all three types of neurons (Fig. S10D-G).
To further reduce contributions of spurious, random low-level expression, we only consider those receptors with at least 10 raw counts per cell and are found in at least 5 cells to evaluate co-expression among receptors. We split cells into groups according to cell-type, age, and neuronal lineage and plot the percentage of cells expressing zero, one, two, and three or more species (Fig. 5B-C). This analysis shows that receptor expression specificity increases as the lineages progressed from progenitors into mature neurons. Immature neurons have more cells co-expressing receptors than mature cells. More co-expressions are observed in the younger animals than the older ones. This line of evidence indicates that the co-expression we have found is not from experimental artifacts but reflects real biological events. Contamination would not be selectively enriched in immature cells and not present in the INP cells.
We performed Fisher’s exact test using contingency tables for every pairing of expressed receptor genes. For the pairs that pass the test, we generated Circos plots to show the genomic loci for all significant V1R, V2R, and interclass pairs (Fig. 5D-F). This analysis showed that 47.7% of co-expressed receptors are co-localized on the same chromosome.
We found a few sets of co-expressed VRs that would have strong implications for how pheromone signals are detected. Vmn1r85 (V1rj3) and Vmn1r86 are two receptors located next to each other on Chr. 7, sharing a high level of homology, and belonging to the V1rj clade. They have a high level of co-expression but are not co-expressed with Vmn1r89 (V1rj2), located ∼100Kb away, even though both Vmn1r85 and Vmn1r89 receptors are activated by sulfated estrogen and carry information about the estrus status of mature female mice63,67,100 (Fig. 5D). Another set of receptors that recognize female-specific pheromone cues are the V1re clade receptors. We found that Vmn1r185 (V1re12), which recognizes female identity pheromones, was co-expressed with its close homolog Vmn1r184 gene, which is about 350kb away on Chr 7. They are not co-expressed with Vmn1r69 (V1re9), which also recognizes female pheromones, but is located 16Mb apart on Chr. 763,100,116,117 (Fig. 5D).
On Chr. 7 there are two other major clusters of V1R genes that show co-expression. One cluster includes Vmn1r55-Vmn1r64, 10 genes belonging to the V1rd clade and located within a 600Kb region. Another one includes Vmn1r167, Vmn1r168, and Vmn1r169, which appear to be specifically paired with Vmn1r175, Vmn1r177, and Vmn1r176, respectively. These receptor pairs are arranged in a 300Kb region with a head-head orientation. Outside of Chr. 7, several small clusters on Chr. 6 and one large cluster on Chr. 3 exhibit significant co-expression of V1Rs.
The V2R neurons coordinately express one common V2R and one specific receptor118,119. Our co-expression analysis confirms the broad association of Vmn2r1-7, which are located on Chr. 3, with other receptors across various chromosomal locations (Fig. 5E). In addition, we also identified co-expression patterns among the specifically expressed V2Rs. Among these receptor genes, intra-chromosomal co-expressions are observed for receptors residing on Chr. 7, Chr. 17, and Chr. 5. There is also inter-chromosomal co-expression between one locus on Chr. 17 with a cluster on Chr. 7 (Fig. 5E).
Lastly, we have observed co-expression of receptors across different classes of receptors (Fig. 5F). Notably, Vmn2r1-3 are co-expressed with several Vmn1r receptors on Chr. 3. Vmn2r7 is co-expressed with a group of Vmn1r genes on Chr. 13. The odorant receptor Olfr344 is co-expressed with several V1R and V2R genes.
In past studies, VR expression was examined using in situ hybridization, immunostaining, or genetic labeling. The traditional histological methods were not sensitive enough to quantitatively measure signal strength. Moreover, pairwise double in situ is too laborious to capture co-expression of two or more receptors. To verify co-expression of VR genes by individual VSNs, we selected 30 VR genes based on scRNA-seq data and used Molecular Cartography© to examine their expression patterns in situ.
Given the high incidents of co-expression between receptors in theVmn1r55-64 cluster on Chr. 7 (Fig. 5D), we included 5 probes for this set of genes and confirmed colocalization of pairs in the VSNs (Fig. 5G), Interestingly, we did not find cells expressing more than two receptor genes for this set. We also confirmed the co-expression between Vmn1r85 and Vmn1r86 as indicated by single cell data (Fig. 5D and G).
We confirmed the co-expression among genes from the V2R genes (Fig. 5H)71,120. Consistent with previous reports71,72, Vmn2r2, Vmn2r3, Vmn2r6, and Vmn2r7 are comingled in several cells, but Vmn2r1 is not co-expressed with these 4 broadly expressed V2Rs (Fig. 5H). Outside of the Vmn2r1-7 group, we find that Vmn2r81 is co-expressed with Vmn2r20 or Vmn2r24 but without any of the Vmn2r1-7 transcripts (Fig. 5H). We detected more incidents of co-expression between Fpr3 and Fpr-rs4, and between Fpr3, Fpr-rs3, Fpr-rs4 with V1Rs than with V2Rs (Fig. 5I and J). We also found co-expressions that are not predicted by the single cell analysis (Fig. 5K). The discrepancy likely can be attributed to the relatively low-level expression of one of the receptor genes. This type of co-expression may not pass the strict criteria we set for the single cell analysis.
A surface molecule code for individual receptor types
VSNs expressing a given receptor type project to the AOB to innervate glomeruli distributed in quasi-stereotypical positions69,73,103. The high number of glomeruli innervated by a given VSN type raises the question about mechanisms that specify the projection patterns and the connection between the VSNs and the mitral cells. For a genetically specified circuit that transmits pheromone and other information to trigger innate behavioral and endocrine responses, there must be molecules that instruct specific connections among neurons. Several studies have revealed the requirement of Kirrel2, Kirrel3, Neuropilin2 (Nrp2), EphA, and Robo/Slit in vomeronasal axon targeting to the AOB121–127. However, how individual guidance molecules or their combinations specify connectivity of individual VSN types is completely unknown. Here we leverage the unbiased dataset to identify surface molecules that may serve as code for circuit specification.
We identified 307 putative axon guidance (AG) molecules, including known cell surface molecules involved in transcellular interactions and some involved in modulating axon growth. Using this panel, we calculated pairwise similarity between two VR genes, and the similarity in their guidance molecule expression. Analysis of the relationship indicates the there is a general increase in guidance molecule similarity associated with VR similarity (Fig. 6A). Consistent with this observation, when we plotted the similarity of surface molecule expression among cells expressing different receptors, we found islands of similarity among the receptor pairs (Fig. 6B). We then conducted correlation analysis between the guidance molecules with the V1Rs (Fig. 6C) and V2Rs (Fig. 6D). Consistent with previous reporting, we found that Kirrel2 was associated with nearly half of the V1Rs and Kirrel3 was mainly associated with V2Rs that project to the caudal AOB121,122. Robo2 is associated with nearly all V2Rs. We found that Teneurin (Tenm2 and Tenm4) and protocadherin (Pcdh9, Pcdh10, and Pcdh17) genes were associated with specific receptors128,129. Epha5, Pdch10, Tenm2, Nrp2 are also strongly associated with V1Rs with partial overlap with each other. Pcdh9, Tenm4, Cntn4, EphrinA3, Pchd17, as well as Kirrel2 and Kirrel3 all show association with specific V2Rs. Numerous guidance molecules that are not broadly expressed are also associated with individual VRs.
Based on these associations, and supported by existing literature, we propose a model of the specification of separate groups of VRs. In this model, Robo2 separates the rostral vs. caudal AOB. Robo2+ V2R neurons project to the caudal AOB whereas Robo2- V1R cells project to the rostral AOB126,127,130. Among V1Rs, Kirrel2 expression distinguishes between two groups of cells121,131. The Kirrel2+ population can be further separated into EphA5+ and EphA5-population122,131. The EphA5+ population can be separated further by Pcdh10 expression. In the Kirrel2-cells, EphA5, Pcdh10, Tenm4, and Tenm2 mark separate groups. Ncam1, EphA5, and Cntn4 may contribute to specifying small sets of cells. For the V2Rs, Pcdh9, Cntn4, Tenm4, Tenm2, and Pcdh17 have a decreasing range of expression, which may be used to specify increasingly refined connection. Notably, even though Kirrel2 and Pcdh10 are mostly detected in the V1Rs, they are also expressed by small sets of the V2Rs.
Transcriptional regulation of receptor and axon guidance cues
The specification of a vomeronasal circuit needs to be tied to receptor expression. We next address whether a transcriptional code is associated with both VR and AG molecule expression. We calculated pairwise similarity among the receptors according to their expression of TF and AG genes (Fig. 7A). The result shows that V1R and V2R are distinctively separated. Besides a small group of broadly expressed V2Rs, all VR types are unique in their gene expression. Fpr types are more similar in their expression profile with the V1Rs, but the OR types are intermingled with both VR types.
To further determine the TF/AG associations that are specific to individual receptor types, we applied stringent criteria to calculate the Jaccard Index for each pair of TF/AG (Fig. 7B), which reflect the statistical probability of co-expression within the same cell. The analysis reveals unique combinations for every receptor type (Fig. 7C-G). Even though the only TF associated with V1R expression is Meis2, other TFs are involved in specifying AG gene expression. For example, cells expressing V1Rs with high sequence homology and in chromosomal proximity share a similar set of TF and AG genes, but the expression patterns are distinct from each other (Fig. 7D). For broadly expressed V2Rs, Vmn2r3 and Vmn2r7, which are co-expressed by the same set of cells, share nearly identical TF/AGs (Fig. 7E). In contrast, Vmn2r1, which does not co-express with either Vmn2r3 or Vmn2r7, lacks the expression of Pou2f1 and Tenm2 despite sharing all other genes. Other V2R types also show distinct TF/AG associations (Fig. 7F-G).
Discussion
Single cell RNA-seq analysis provides an unprecedented opportunity to identify cell types and determine genes associated with individual cells, but it is not without pitfalls. At the current state of the art, the depth of sequencing only allows sampling of transcripts that are expressed at relatively high levels. Sequencing dropouts and potential contamination also can complicate the analysis. Using the current state of the art tools, and applying conservative criteria, we provide an in-depth look at the molecular and cellular organization of the mouse VNO. The analyses reveal new cell types, specific co-expression of receptors, and transcriptional regulation of lineage specification. Moreover, our analyses uncover specific associations between transcription factors, surface guidance molecules, and individual receptor types that may determine the wiring specificity in the vomeronasal circuitry.
The molecular distinction between the sVSNs and the classic VSNs indicates that they may serve a specific function. They are different from the solitary chemosensory cells that are trigeminal in nature 132. We speculate that these cells may secrete olfactory binding proteins and mucins in response to VNO activation. The VNO is a semi-blind tubular structure. Chemical cues are actively transported into the VNO and can only be cleared by active transport systems, which are thought to be carried out by the lipocalin family of proteins, or by degradation 132–135. These proteins are important to protect the integrity of neuroepithelia. They are generally produced by the sustentacular cells or the Bowman’s gland 135. It is plausible that the sVSNs can produce specific lipocalin proteins based on the ligand detected by the VRs they express. These cells may also convey chemosensory information to the AOB, but we do not know if they project to the central brain. We also have identified a class of canonical OSNs in the VNO. Previous reports show that these neurons project to AOB. It is plausible that the OSNs can detect a set of volatile odors that carry species-specific information and directly convey it to brain areas that regulate innate responses. Our list of these ORs could direct effort to identify these odors to reveal their ethological relevance.
The co-expression of multiple VRs in individual VSNs is intriguing, as a previous analysis of the MOE detected minimal co-expression of ORs21. Importantly, there is a higher propensity for VR co-expression between certain receptor pairs. Notably, there is co-expression of receptors sharing common ligands, or that are similar in their sequences. These observations indicate that co-expressed receptors may serve redundant function detecting the cues. For example, the V1rj receptors are cognate receptors for sulfates estrogen and estrus signals. Their co-expression indicates that the neurons detect the same class of molecules redundantly. Moreover, co-expression of similarly tuned receptors makes it plausible for heterotypic convergence, when these neurons converge into glomeruli that express one or the other receptors.
How specificity of neuronal connections in the olfactory system is determined remains unknown. In the main olfactory system, glomerular positions are coarsely specified along the anterior-posterior and dorsal-ventral axes by gradients of axon guidance molecules whereas the sorting of axons according to the odorant receptors is mediated by homophilic attraction and heterotypic repulsion using a different set of guidance molecules 136. Spontaneous neural activities determine the expression of both sets of molecules. It is not known whether VSNs utilize the same mechanism. Given the multi-glomerular innervation patterns by the VSNs, it is exceedingly difficult to determine the contribution of individual guidance molecules to specifying VSN innervation.
We have identified guidance molecules associated with individual VRs that potentially constitute a code set that specifies VSN axon projections and their connection with postsynaptic cells. Each receptor type has a unique combination of guidance molecules expressed, which provides a basis for axon segregation and convergence. There are a few molecules that are shared broadly by various VSN types. These can be used to instruct general spatial locations of the VSN axons. For example, Robo2 separates the anterior vs. posterior AOB. Knockout of Robo2 causes mistargeting of V2R neurons to the rostral AOB. Our models also indicate that Kirrel2 and Kirrel3 are expressed by nearly half of the VR types in partially overlapping patterns. Deletion of Kirrel2 or Kirrel3 leads to disorganization of glomeruli in the posterior AOB. Protocadherins and tenurins add new dimensions to this code. We also identified several guidance molecules that are more specifically associated with individual VRs. They could provide additional cues to separate axons that share broadly expressed guidance molecules.
We have identified lineage relationships among cells in the VNO and a dynamic transcriptional cascade that likely specifies cell types during development. While our model agrees with that of Katreddi et al26 on the main transcription factors that specify the lineage, it adds more details on both the induction and suppression of genes in specifying the cell fate. For example, we confirm Meis2 and Tfap2e as transcription factors that maintain the V1R and V2R fate, but we also found that the down regulation of Neurog1, but not NeuroD1, is associated with a transition from early INPs to late INPs. The downregulation of Sp8, Nfib, and Bcl11b is likely important for committing to the V1R lineage for late INPs. On the other hand, downregulation of Sp8 and upregulation of Fezf1, Olig2, and Tshz2 likely set up commitment to the OSN fate. We also find that in all three lineages, the expression of Tshz2 is associated with transition to the immature neuronal fate from the late INPs.
We observed a striking difference between V1R and V2R VSNs in the transcription factors associated with receptor choice. There is no overt association between V1R with specific transcription factors. This observation is reminiscent of OSNs in the olfactory epithelium, where OR expression is stochastic and mediated by de-repression of epigenetically silenced OR loci106–111. The absence of specific transcription factors with individual V1R choice suggests that a similar mechanism may operate in the V1R VSNs. Monoallelic expression of V1Rb2 supports this notion69. On the other hand, we observed that for individual V1R types, there are specific associations between transcription factors with guidance molecules. This observation implies that the expression of guidance molecules is determined by combinations of transcription factors even though these transcription factors may not determine V1R expression. This is also reminiscent of the OSNs, where the expression of guidance molecules is determined by spontaneous neural activities137–139. That is, once the receptor choice is made, the specific receptor being expressed determines the guidance molecules to specify their projection patterns.
In direct contrast, V2R VSNs likely use combinations of transcription factors to specify receptor expression as well as guidance molecules. Some of the transcription factors that we observe to be associated with V2R expression are also associated with guidance molecule expression. For example, Pou2f1, Atf5, and Zfp268 are involved in both processes.
Acknowledgements
We thank McKenzie Treese, KyeongMin Bae, Fang Liu, and members of Lab Animal Support Facility at Stowers for their technical support. We would like to acknowledge the University of Kansas Medical Center’s Genomics Core for their support in generating data on the Illumina NovaSeq 6000 System. The core is supported by the following grants: Kansas Intellectual and Developmental Disabilities Research Center (NIH U54 HD 090216), the Molecular Regulation of Cell Development and Differentiation – COBRE (P30 GM122731-03) and the NIH S10 High-End Instrumentation Grant (NIH S10OD021743). This work was supported by fundings from NIH (R01DC008003 and R01DC020368) and Stowers Institute for Medical Research to C.R.Y.
Additional information
Author contributions
L.M. and C.R.Y. conceived the study. A.F. and L.M. performed single cell RNA-seq experiments. T.C. and S.M. prepared samples for Molecular Cartography analysis. A.P. and A.S. facilitated single cell library preparation and RNA sequencing. M.H. performed bioinformatic analysis. M.H and T.C created figures and tables. C.R.Y. acquired funding, supervised the project, and wrote the final manuscript.
Declaration of interests
The authors declare no competing interests.
STAR Methods
Key resources table
Resource availability
Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by C. Ron Yu (ryu@stowers.org).
Data and code availability
All RNA-seq data are available from the NCBI GEO server (GSE252365). All original data generated in this study will be available for download at Stowers original data repository upon publication. No custom generated computer code was used for analysis.
An HTML file containing relevant figures and statistics from the study, as well as tables showing co-expression and differential expression results, can be accessed in the following URL: https://ronyulab.github.io/VNO-Atlas/.
Experimental model and subject details
Wildtype CD1 postnatal day 14 (P14) pups and adults (P56) were used for the experiment. Both sexes were randomly assigned to the experiment. All animals were maintained in Stowers LASF with a 14:10 light cycle and provided with food and water ad libitum. Experimental protocols were approved by the Institutional Animal Care and Use Committee at Stowers Institute and in compliance with the NIH Guide for Care and Use of Animals.
Methods
scRNA library preparation & sequencing
Mice VNOs were dissected in cold oxygenated ACSF following Ma et al, 2011140. Dissected epitheliums were dissociated in papain solution (20mg/mL papain and 3mg/mL L-cysteine in HBSS) with RNase-free DNase I (10unit) at 37°C for 15-20 minutes. 0.01% BSA in PBS was added to the digestion solution before filtering with 70µm and 30µm filters (pluriSelect). Dissociated cells were washed twice in 0.01% BSA with final volume 1mL, followed by Draq5 (25µM) and DAPI (0.5µg/mL) staining 5min on ice. Draq5+/DAPI- cells (live/nucleated cells) were sorted on BD Influx cytometer (BD Bioscience) with 100µm nozzle. Dissociated sorted cells were assessed for concentration and viability via Luna-FL cell counter (Logos Biosystems). Cells deemed to be at least 90% viable were loaded on a Chromium Single Cell Controller (10x Genomics), based on live cell concentration. Libraries were prepared using the Chromium Next GEM Single Cell 3’ Reagent Kits v3.1 (10X Genomics) according to manufacturer’s directions. Resulting cDNA and short fragment libraries were checked for quality and quantity using a 2100 Bioanalyzer (Agilent Technologies) and Qubit Fluorometer (Thermo Fisher Scientific). With cells captured estimated at ∼5,500-8,000 cells per sample, libraries were pooled and sequenced to a depth necessary to achieve at least 40,000 mean reads per cell on an Illumina NovaSeq 6000 instrument utilizing RTA and instrument software versions current at the time of processing with the following paired read lengths: 28*8*91bp.
scRNA-Seq pre-processing and QC-filtering
Gene-by-cell barcode count matrices were generated from raw FASTQ files using the kallisto|bustools (v0.48.0|v0.41.0)141,142 workflow with the Mus musculus genome assembly GRCm39 (mm39) and GTF gene annotation files retrieved from ENSEMBL release 104143. All downstream QC-filtering and analysis was performed in an R environment (v4.0.3)144. Empty droplets were estimated and filtered from the data using the barcodeRanks function of the DropletUtils package (v1.10.3)145 with the lower bound of UMI-counts set to 100. The count data was then imported into R (v4.0.3) using the Seurat package (v4.3.0)146 and ambient RNA contamination was automatically estimated with the autoEstCont function and removed with the adjustCounts function from the SoupX package (v1.5.2)115. Cell barcodes representing multiplets were identified and removed with Scrublet (v0.2.3)147 interfacing with Python using reticulate (v1.30)148. Cell barcodes expressing <750 genes, >2.5 standard deviations above the mean number of genes or counts, or >5% of reads originating from mitochondrial genes were removed from the downstream analysis.
scRNA-Seq integration and clustering
To cluster cells from multiple samples, raw gene counts for each sample were normalized with the SCTransform function (v0.3.2)149 using the glmGamPoi method from the glmGamPoi package (v1.6.0)150, samples were then integrated using the Seurat integration pipeline. After principal component analysis (PCA) the dimensionality of the whole integrated VNO dataset was estimated to be 25 based on visual identification of an ‘elbow’ using the output from the ElbowPlot function. The shared nearest neighbor graph was constructed with the FindNeighbors function. Then, using the FindClusters function, an optimal resolution of 0.7 was chosen for the discovery of broad cell types in the VNO by running the clustree function from the clustree package (v0.30)151 to evaluate cluster-level mean expression for a curated list of VNO cell-type marker-genes in a range of resolutions between 0 and 1, with intervals of 0.1; optimal resolution was considered the lowest resolution at which the expression of Neurod1 and Ascl1 diverge into two clusters, as these genes are unique marker-genes for immediate neural progenitors (INPs) and globose basal cells (GBCs), respectively. At 0.7 resolution 31 clusters were observed.
Differential gene expression analysis and cell labeling
Clusters were assigned to cell-types using a list of marker-genes: Omp was used to broadly mark mature neurons, Gnai2, Gng13, and Meis2 marked the V1R subtype, Gnao1, Robo2, and Tfap2e marked the V2R subtype, Gnal and Cnga2 marked OSNs, Gap43, Stmn2, Bcl11b, and Lhx2 marked immature neurons, Neurod1 and Neurog1 marked INPs, Ascl1 and Ccnd1 marked GBCs, Krt5 and Krt14 marked HBCs, Sox9 and Hepacam2 marked MVs, Fezf2 and Sox2 marked Sustentacular cells, S100b and Plp1 marked OECs, Cd34 and Cdh5 marked Endothelial cells, Acta2, Col1a2, and Mgp marked Lamina Propria cells, Dock2 marked immune cells, Cx3cr1 and Ctss marked Microglia, and Il7r and Trbc2 marked T-cells. First, we ran the FindAllMarkers function on the clustered normalized counts, limited to genes present in ≥ 50% of each cluster’s cells with an absolute value log2 fold-change ≥ 0.5; additionally, we ran the FeaturePlot function to plot expression of the marker-genes in 2D UMAP space. Clusters were then manually assigned to one cell-type both by visually ascertaining marker-gene expression overlap with cluster identity and by statistically validating significant marker-gene enrichment within cluster.
Neuronal Lineage
Cells labelled as GBC, INP, iVSN, iOSN, mVSN, mOSN, or sVSN were subset from the integrated whole VNO dataset, split by animal sample, then re-integrated and re-clustered, as described above. The FindNeighbors function was run using 12 PCs and the FindClusters function was run with a resolution of 6.0 which resulted in 84 clusters. Clusters were assigned to cell-types as described previously, while a further distinction between early and late INPs was inferred from the expression of Ascl1, Neurod1, Neurog1, and Gap43 in 2D UMAP space. Differential gene expression analysis was performed between the novel sVSN cluster and the mature V1R, V2R, and OSN clusters, independently, running the FindMarkers function with the logfc.threshold and min.pct parameters set to 0 to accommodate downstream gene set enrichment analysis (GSEA).
Gene set enrichment analysis
Gene ontology (GO) terms from the biological process, molecular function, and cellular compartment categories along with their associated gene sets were retrieved with the msigdbr function from the msigdbr package (v7.5.1)152. Ranked Wald test differential expression results between sVSNs and V1R/V2R mVSNs, respectively, were input into the stats parameter of the fgsea function from the fgsea package (v1.20.0)153,154 along with the GO term gene sets. The top significant GO terms were plotted using the ggplot function from the ggplot2 package.
Sex and age differences
To examine broad differences in gene expression between male and female mice and between P14 and P56 mice, we ran the FindMarkers function on the normalized count data using all cells in the neuronal lineage, all genes present in the data, and no threshold on log2 fold-change. Significant differential gene expression results (padj ≤ 0.05) from the male/female and the P14/P56 test were used for GSEA, and the results were plotted using the ggplot function from the ggplot2 package.
Immediate neural progenitors and immature vomeronasal sensory neurons
Cells previously identified as early and late INP, iVSN, iOSN, or mOSN, were subset from the neuronal dataset, split and reintegrated, as above, using 15 PCs and a resolution of 2.5, resulting in 30 clusters. Differential gene expression analysis was performed with the FindMarkers function between two clusters showing either V1R or V2R like properties but previously identified uniformly as early INPs in the whole neuronal lineage dataset. The FeaturePlot function, from the Seurat package, was used to show normalized expressions for genes of interest.
Trajectory inference and differential gene expression analysis
V1R/V2R lineage determination
To explore transcriptional differences over pseudotime between V1R and V2R VSNs, we subset GBCs, early and late INPs, V1R and V2R iVSNs, and mVSNs from the neuronal dataset, split the data by sample and reintegrated, then performed PCA. Trajectory inference analysis was performed with the slingshot function from the slingshot package (v1.8.0)99 using the first 12 PCs and cell type labels from the neuronal dataset, with an input starting cluster of GBCs and two input end clusters for V1R and V2R mVSNs. To determine the nknots parameter for the fitGAM function from the tradeSeq package (v1.4.0)155, we ran the evaluateK function with the raw count matrix, and the pseudotime and cell weight values output from slingshot. We then ran the fitGAM with nknots=5. We then ran the patternTest function to test for differences in gene expression patterns over pseudotime between the V1R and V2R lineages.
Pseudotime analysis of sVSNs
To determine whether sVSNs represent an immature version of canonical VSNs, we performed trajectory inference analysis on the neuronal lineage with only cells belonging to the OSN lineage removed. Using the slingshot function we input the first 5 PCs and set the starting cluster to GBCs and three end-clusters to V1R mVSNs, V2R mVSNs, and sVSNs, respectively. Pseudotime values were assigned to cells based on their lineage membership and were subsequently plotted with the FeaturePlot function.
Gene co-expression analysis
VR co-expression
To investigate cell-level diversity of vomeronasal receptor species across all cells in the neuronal dataset we calculated the Shannon diversity index for the raw gene counts for all Vmn1r, Vmn2r, Olfr, and Fpr genes using the diversity function from the vegan R package (v2.6-4)156 with default parameters. To determine what proportion of cells in the neuronal lineage had one, two, or three or more receptor species, we set a threshold of ≥10 raw counts for a receptor to be considered ‘present’. To test whether co-expressing vomeronasal receptor species were significant, using the same raw-count threshold of ≥10, we gathered a list of all cell barcodes where a receptor was observed, for all receptors. Using the lists of cell barcodes associated with the receptors, we ran the newGOM function from the GeneOverlap R package (v1.26.0)157, which calculates p-values using Fisher’s exact test on a contingency table. P-values were then corrected for multiple testing using the Benjamini-Hochberg procedure.
Using the circlize R package (v0.4.15)158, we plotted all co-expressed receptor pairs on circos plots showing each receptors genetic location and the number of cells expressing the pair, for all significant pairings (padj ≤ 0.05).
VR co-expression with axon guidance (AG) and transcription factor (TF) genes
To ascertain vomeronasal receptor co-expression with AG genes expressed at the plasma membrane, and with DNA binding TF genes, respectively, we used the Mouse Genome Informatics database to find all genes associated with the biological process gene ontology (GO) term ‘axon guidance’, or with the molecular function GO term ‘DNA-binding transcription factor activity’. For the axon guidance gene set, we subset genes that were expressed at the plasma membrane. We set the VR raw count threshold to ≥10 counts and the AG and TF gene raw count threshold to ≥3. Then, using the cell barcodes associated with each VR and the cell barcodes associated with candidate AG and TF genes, we ran the newGOM function to find significant co-expression (padj ≤ 0.05) for all VRxAG and VRxTF pairs.
VR-specific co-expression of AG and TF genes
To test whether AGs and TFs co-expressed for a given VR, we gathered all cell-barcodes where there was significant VRxAG or VRxTF co-expression. Then, using the same contingency table scheme as above, we looked for significant co-expression (padj ≤ 0.05) between all AGs and TFs previously found to co-express with a given VR.
Spatial transcriptomics
Samples
VNO tissues were dissected from 7-8 weeks old C57BL/6J mice. Briefly, mice were anesthetized with urethane at a dose of 2000 mg/kg body weight. Following general anesthesia, mice heads were decapitated, and the lower jaw was removed by cutting the mandible bone with scissors. The ridged upper palate tissue was peeled off to expose the nasal cavity. A surgical blade was inserted between the two upper incisors to expose the VNO. The whole VNO was carefully extracted by holding onto the tail bone and slowly lifting it up from the nasal cavity. The dissected VNO was immediately transferred to cold 1X PBS on ice, and subsequently embedded in frozen section media (Leica Surgipath FSC 22, Ref # 3801481) and frozen on liquid nitrogen. Frozen samples were sectioned at 10 µm thickness using the Thermo Scientific CryoStar NX70 cryostat. VNO sections were placed within capture area of cold slides that were provided by Resolve Biosciences. Slides were sent to Resolve Biosciences on dry ice for spatial transcriptomics analysis. Resolve Molecular Cartography© protocols remain proprietary and were not disclosed. The probe design, tissue processing, imaging, spot segmentation, and image preprocessing were all performed using the Resolve Biosciences platform. Names and ENSEMBL IDs for genes probed are available in this study’s public repository at https://ronyulab.github.io/VNO-Atlas/.
Analysis
Regions of interest in the VNO were selected on brightfield images provided by Resolve Biosciences. Cell segmentation on final images was performed in QuPath 159. Detected gene transcripts were then assigned to the segmented cells, thereby creating a gene-count matrix for each sample. To predict cell types, count matrices for all samples were imported into R then normalized using the SCTransform function with the glmGamPoi method. The samples were then integrated using the Seurat integration pipeline. Both the integrated whole VNO scRNA-seq dataset and the integrated Resolve molecular cartography dataset were renormalized with SCTransform with the default method using ncells=3000, then RunPCA was called on the renormalized data. Using the whole VNO scRNA-seq dataset as a reference and the molecular cartography dataset as a query, we ran the FindTransferAnchors function, then we ran the TransferData function to create a table of prediction score values for each cell in the spatial dataset. Cells were then labeled by type using the maximum prediction score for each cell. Images showing co-localization of vomeronasal receptors were obtained in ImageJ using genexyz Polylux (v1.9.0) tool plugin from Resolve Biosciences.
Region of neurogenesis
To test the hypothesis that neurogenesis occurs in the marginal zone of the VNO we first set a minimum cell-type prediction-score threshold of ≥ 0.3; all cells below the threshold were labeled “unknown”. Then we used the simple features R package, sf (v1.0.16) to delineate regions of interest in the VNO. We excluded the non-neuronal region of the VNO from the analysis. Using the intersectional boundary between neural and non-neuronal epithelia as the center, we quantified the cells falling within a 750-pixel radius as in the marginal zones.
Those fall out of the 750-pixel but within a 1500-pixel radius were quantified as in the intermediate zones. All remaining cells were classified as occurring in the “main zone”. Plots displaying the results were created using the function, ggplot.
References
- 1.PheromonesAmerican Elsevier Pub. Co.
- 2.Pheromones and animal behaviour : communication by smell and tasteCambridge University Press
- 3.Pheromonal communication in vertebratesNature 444:308–315https://doi.org/10.1038/nature05404
- 4.Molecular organization of vomeronasal chemoreceptionNature 478:241–245https://doi.org/10.1038/nature10437
- 5.Encoding gender and individual information in the mouse vomeronasal organScience 320:535–538https://doi.org/10.1126/science.1154476
- 6.Pheromones and behavior in miceActa Neurol Psychiatr Belg 69:529–538
- 7.Acceleration and delay of puberty in female housemice: methods of delivery of the urinary stimulusDev Psychobiol 14:487–497
- 8.Structure and function of the vomeronasal system: an updateProg Neurobiol 70:245–318
- 9.Pheromones and reproduction in mammalsAcademic Press
- 10.Coordination of social signals and ovarian function during sexual developmentJ Anim Sci 67:1841–1847
- 11.Sexual behavior and aggression in male mice: involvement of the vomeronasal systemThe Journal of neuroscience : the official journal of the Society for Neuroscience 4:2222–2229
- 12.Mediation of male mouse urine marking and aggression by the vomeronasal organPhysiology & behavior 37:655–657
- 13.Vomeronasal functionChemical senses 23:463–466
- 14.Sensory, hormonal, and neural control of maternal aggression in laboratory rodentsNeurosci Biobehav Rev 26:869–888
- 15.The neuroendocrine basis of social recognitionFront Neuroendocrinol 23:200–224https://doi.org/10.1006/frne.2002.0229
- 16.Cyclic regulation of sensory perception by a female hormone alters behaviorCell 161:1334–1344
- 17.Developmental study on the vomeronasal organ in the rat fetusJournal of Reproduction and Development 39:47–54
- 18.Prenatal development of the mammalian vomeronasal organMicroscopy research and technique 41:456–470
- 19.The human olfactory transcriptomeBMC genomics 17:1–18
- 20.A transcriptional rheostat couples past activity to future sensory responsesCell 184:6326–6343
- 21.Single-cell transcriptomics reveals receptor transformations during olfactory neurogenesisScience 350:1251–1255
- 22.Deconstructing olfactory stem cell trajectories at single-cell resolutionCell stem cell 20:817–830
- 23.A Population of Navigator Neurons Is Essential for Olfactory Map Formation during the Critical PeriodNeuron 100:1066–1082https://doi.org/10.1016/j.neuron.2018.09.051
- 24.Analysis of the vomeronasal organ transcriptome reveals variable gene expression depending on age and function in rabbitsGenomics 113:2240–2252
- 25.Pronounced strain-specific chemosensory receptor gene expression in the mouse vomeronasal organBMC Genomics 18https://doi.org/10.1186/s12864-017-4364-4
- 26.Notch signaling determines cell-fate specification of the two main types of vomeronasal neurons of rodentsDevelopment 149
- 27.Sociosexual behavior requires both activating and repressive roles of Tfap2e/AP-2epsilon in vomeronasal sensory neuronseLife 11https://doi.org/10.7554/eLife.77259
- 28.A novel family of putative pheromone receptors in mammals with a topographically organized and sexually dimorphic distributionCell 90:763–773
- 29.A new multigene family of putative pheromone receptorsNeuron 19:371–379
- 30.Formyl peptide receptor-like proteins are a novel family of vomeronasal chemosensorsNature 459:574–577https://doi.org/10.1038/nature08029
- 31.Formyl peptide receptors are candidate chemosensory receptors in the vomeronasal organProc Natl Acad Sci U S A 106:9842–9847https://doi.org/10.1073/pnas.0904464106
- 32.A novel family of genes encoding putative pheromone receptors in mammalsCell 83:195–206
- 33.A multigene family encoding a diverse array of putative pheromone receptors in mammalsCell 90:775–784
- 34.Multiple new and isolated families within the mouse superfamily of V1r vomeronasal receptorsNat Neurosci 5:134–140
- 35.Origin and evolution of the vertebrate vomeronasal system viewed through system-specific genesBioessays 28:709–718https://doi.org/10.1002/bies.20432
- 36.Comparative genomic analysis identifies an evolutionary shift of vomeronasal receptor gene repertoires in the vertebrate transition from water to landGenome Res 17:166–174https://doi.org/10.1101/gr.6040007
- 37.Vomeronasal receptors in vertebrates and the evolution of pheromone detectionAnnual review of animal biosciences 5:353–370
- 38.Species specificity in rodent pheromone receptor repertoiresGenome Res 14:603–608https://doi.org/10.1101/gr.2117004
- 39.Dynamic evolution of V1R putative pheromone receptors between Mus musculus and Mus spretusBMC Genomics 10https://doi.org/10.1186/1471-2164-10-74
- 40.Odorant and vomeronasal receptor genes in two mouse genome assembliesGenomics 83:802–811
- 41.Electrophysiological characterization of chemosensory neurons from the mouse vomeronasal organThe Journal of neuroscience : the official journal of the Society for Neuroscience 16:4625–4637https://doi.org/10.1523/JNEUROSCI.16-15-04625.1996
- 42.TRICK or TRP? What Trpc2(-/-) mice tell us about vomeronasal organ mediated innate behaviorsFrontiers in neuroscience 9https://doi.org/10.3389/fnins.2015.00221
- 43.Evidence for different chemosensory signal transduction pathways in olfactory and vomeronasal neuronsBiochem Biophys Res Commun 220:900–904https://doi.org/10.1006/bbrc.1996.0503
- 44.Hyperpolarization-activated cyclic nucleotide-gated channels in mouse vomeronasal sensory neuronsJ Neurophysiol 100:576–586https://doi.org/10.1152/jn.90263.2008
- 45.Calcium-activated chloride channels in the apical region of mouse vomeronasal sensory neuronsThe Journal of general physiology 140:3–15https://doi.org/10.1085/jgp.201210780
- 46.Conditional knockout of TMEM16A/anoctamin1 abolishes the calcium-activated chloride current in mouse vomeronasal sensory neuronsThe Journal of general physiology 145:285–301https://doi.org/10.1085/jgp.201411348
- 47.Odors activate dual pathways, a TRPC2 and a AA-dependent pathway, in mouse vomeronasal neuronsAm J Physiol Cell Physiol 298:C1253–1264https://doi.org/10.1152/ajpcell.00271.2009
- 48.Calcium-activated chloride current amplifies the response to urine in mouse vomeronasal sensory neuronsThe Journal of general physiology 135:3–13https://doi.org/10.1085/jgp.200910265
- 49.Arachidonic acid plays a role in rat vomeronasal signal transductionThe Journal of neuroscience : the official journal of the Society for Neuroscience 22:8429–8437
- 50.A diacylglycerol-gated cation channel in vomeronasal neuron dendrites is impaired in TRPC2 mutant mice: mechanism of pheromone transductionNeuron 40:551–561
- 51.Pheromonal recognition memory induced by TRPC2-independent vomeronasal sensingEur J Neurosci 23:3385–3390https://doi.org/10.1111/j.1460-9568.2006.04866.x
- 52.Sensory transduction in vomeronasal neurons: evidence for G alpha o, G alpha i2, and adenylyl cyclase II as major components of a pheromone signaling cascadeThe Journal of neuroscience : the official journal of the Society for Neuroscience 16:909–918
- 53.Paradoxical contribution of SK3 and GIRK channels to the activation of mouse vomeronasal organNat Neurosci https://doi.org/10.1038/nn.3173
- 54.Requirement of calcium-activated chloride channels in the activation of mouse vomeronasal neuronsNat Commun 2
- 55.Central role of G protein Gαi2 and Gαi2+ vomeronasal neurons in balancing territorial and infant-directed aggression of male miceProceedings of the National Academy of Sciences 116:5135–5143
- 56.G protein G(alpha)o is essential for vomeronasal function and aggressive behavior in miceProc Natl Acad Sci U S A 108:12898–12903https://doi.org/10.1073/pnas.1107770108
- 57.Ultrastructural localization of G-proteins and the channel protein TRP2 to microvilli of rat vomeronasal receptor cellsJ Comp Neurol 438:468–489
- 58.The TRPC2 ion channel and pheromone sensing in the accessory olfactory systemNaunyn Schmiedebergs Arch Pharmacol 371:245–250https://doi.org/10.1007/s00210-005-1028-8
- 59.Loss of sex discrimination and male-male aggression in mice deficient for TRP2Science 295:1493–1500https://doi.org/10.1126/science.1069259
- 60.Altered sexual and social behaviors in trp2 mutant miceProc Natl Acad Sci U S A 99:6376–6381https://doi.org/10.1073/pnas.082127599
- 61.TRPC2 and the Molecular Biology of Pheromone Detection in MammalsTRP Ion Channel Function in Sensory Transduction and Cellular Signaling Cascades Boca Raton: CRC Press/Taylor & Francis
- 62.Ultrasensitive pheromone detection by mammalian vomeronasal neuronsNature 405:792–796https://doi.org/10.1038/35015572
- 63.Tuning properties and dynamic range of type 1 vomeronasal receptorsFrontiers in neuroscience 9https://doi.org/10.3389/fnins.2015.00244
- 64.Distinct signals conveyed by pheromone concentrations to the mouse vomeronasal organThe Journal of neuroscience : the official journal of the Society for Neuroscience 30:7473–7483https://doi.org/10.1523/JNEUROSCI.0825-10.2010
- 65.Pheromone detection mediated by a V1r vomeronasal receptorNat Neurosci 5:1261–1262https://doi.org/10.1038/nn978
- 66.Responses of vomeronasal neurons to natural stimuliScience 289:1569–1572
- 67.Sulfated steroids as natural ligands of mouse pheromone-sensing neuronsThe Journal of neuroscience : the official journal of the Society for Neuroscience 28:6407–6418https://doi.org/10.1523/JNEUROSCI.1425-08.2008
- 68.Gene cluster lock after pheromone receptor gene choiceEMBO J 26:3423–3430https://doi.org/10.1038/sj.emboj.7601782
- 69.Variable patterns of axonal projections of sensory neurons in the mouse vomeronasal systemCell 97:199–208
- 70.Co-expression of putative pheromone receptors in the sensory neurons of the vomeronasal organThe Journal of neuroscience : the official journal of the Society for Neuroscience 21:843–848https://doi.org/10.1523/JNEUROSCI.21-03-00843.2001
- 71.A recent class of chemosensory neurons developed in mouse and ratPLoS One 6
- 72.Combinatorial co-expression of pheromone receptors, V2RsJournal of neurochemistry 103:1753–1763https://doi.org/10.1111/j.1471-4159.2007.04877.x
- 73.A map of pheromone receptor activation in the mammalian brainCell 97:209–220
- 74.Representation and transformation of sensory information in the mouse accessory olfactory systemNat Neurosci 13:723–730https://doi.org/10.1038/nn.2546
- 75.A divergent pattern of sensory axonal projections is rendered convergent by second-order neurons in the accessory olfactory bulbNeuron 35:1057–1066
- 76.A multireceptor genetic approach uncovers an ordered integration of VNO sensory inputs in the accessory olfactory bulbNeuron 50:697–709https://doi.org/10.1016/j.neuron.2006.04.033
- 77.Light microscopic Golgi study of mitral/tufted cells in the accessory olfactory bulb of the adult ratJournal of Comparative Neurology 311:65–83
- 78.Visualizing an olfactory sensory mapCell 87:675–686
- 79.A physicochemical model of odor samplingPLoS Computational Biology 17
- 80.Acquisition of Innate Odor Preference Depends on Spontaneous and Experiential Activities During Critical PeriodbioRxiv
- 81.Molecular Control of Circuit Plasticity and the Permanence of Imprinted Odor MemorybioRxiv https://doi.org/10.1101/2022.03.29.486284
- 82.Regeneration of new neurons is preserved in aged vomeronasal epitheliaJournal of Neuroscience 30:15686–15694
- 83.Reliable sex and strain discrimination in the mouse vomeronasal organ and accessory olfactory bulbThe Journal of neuroscience : the official journal of the Society for Neuroscience 33:13903–13913https://doi.org/10.1523/JNEUROSCI.0037-13.2013
- 84.Multielectrode array recordings of the vomeronasal epitheliumJ Vis Exp https://doi.org/10.3791/1845
- 85.Cells in the vomeronasal organ express odorant receptors but project to the accessory olfactory bulbJournal of Comparative Neurology 498:476–490
- 86.Detection of pup odors by non-canonical adult vomeronasal neurons expressing an odorant receptor gene is influenced by sex and parenting statusBMC biology 14:1–20
- 87.Prenatal development of the mammalian vomeronasal organMicrosc Res Tech 41:456–470https://doi.org/10.1002/(SICI)1097-0029(19980615)41:6<456::AID-JEMT2>3.0.CO;2-L
- 88.Molecular switches in the development and fate specification of vomeronasal neuronsJournal of Neuroscience 31:17761–17763
- 89.Mechanisms underlying pre-and postnatal development of the vomeronasal organCellular and Molecular Life Sciences 78:5069–5082
- 90.Proliferation and migration of receptor neurons in the vomeronasal organ of the adult mouseBrain Res Dev Brain Res 123:33–40
- 91.Neurogenesis in the vomeronasal epithelium of adult rats: evidence for different mechanisms for growth and neuronal turnoverJ Neurobiol 44:423–435
- 92.Smad4-dependent morphogenic signals control the maturation and axonal targeting of basal vomeronasal sensory neurons to the accessory olfactory bulbDevelopment 147
- 93.Bcl11b/Ctip2 controls the differentiation of vomeronasal sensory neurons in miceThe Journal of neuroscience : the official journal of the Society for Neuroscience 31:10159–10173https://doi.org/10.1523/JNEUROSCI.1245-11.2011
- 94.The transcription factor Tfap2e/AP-2ε plays a pivotal role in maintaining the identity of basal vomeronasal sensory neuronsDevelopmental biology 441:67–82
- 95.Co-expression of C/EBPγ and ATF5 in mouse vomeronasal sensory neurons during early postnatal developmentCell and tissue research 378:427–440
- 96.Gli3 regulates vomeronasal neurogenesis, olfactory ensheathing cell formation, and GnRH-1 neuronal migrationJournal of Neuroscience 40:311–326
- 97.Specific mesenchymal/epithelial induction of olfactory receptor, vomeronasal, and gonadotropin-releasing hormone (GnRH) neuronsDev Dyn 239:1723–1738https://doi.org/10.1002/dvdy.22315
- 98.Expression patterns of homeobox genes in the mouse vomeronasal organ at postnatal stagesGene Expression Patterns 21:69–80
- 99.Slingshot: cell lineage and pseudotime inference for single-cell transcriptomicsBMC Genomics 19https://doi.org/10.1186/s12864-018-4772-0
- 100.Integrated action of pheromone signals in promoting courtship behavior in male miceeLife 3https://doi.org/10.7554/eLife.03025
- 101.A single vomeronasal receptor promotes intermale aggression through dedicated hypothalamic neuronsNeuron 110:2455–2469
- 102.A juvenile mouse pheromone inhibits sexual behaviour through the vomeronasal systemNature 502:368–371https://doi.org/10.1038/nature12579
- 103.The male mouse pheromone ESP1 enhances female sexual receptive behaviour through a specific vomeronasal receptorNature 466:118–122https://doi.org/10.1038/nature09142
- 104.Hemoglobin in the blood acts as a chemosensory signal via the mouse vomeronasal systemNature Communications 13
- 105.Allelic inactivation regulates olfactory receptor gene expressionCell 78:823–834
- 106.LHX2-and LDB1-mediated trans interactions regulate olfactory receptor choiceNature 565:448–453
- 107.Co-opting the unfolded protein response to elicit olfactory receptor feedbackCell 155:321–332
- 108.Nuclear aggregation of olfactory receptor genes governs their monogenic expressionCell 151:724–737
- 109.An epigenetic trap stabilizes singular olfactory receptor expressionCell 154:325–336
- 110.Monoallelic expression of olfactory receptorsAnnual review of cell and developmental biology 31:721–740
- 111.Interchromosomal interactions and olfactory receptor choiceCell 126:403–413
- 112.Singular expression of olfactory receptor genesCell 155:274–277
- 113.Gene cluster lock after pheromone receptor gene choiceThe EMBO journal 26:3423–3430
- 114.Clustering of vomeronasal receptor genes is required for transcriptional stability but not for choiceScience Advances 8
- 115.SoupX removes ambient RNA contamination from droplet-based single-cell RNA sequencing dataGigascience 9https://doi.org/10.1093/gigascience/giaa151
- 116.A Molecular Code for Identity in the Vomeronasal SystemCell 163:313–323https://doi.org/10.1016/j.cell.2015.09.012
- 117.Sensory coding mechanisms revealed by optical tagging of physiologically defined neuronal typesScience 366:1384–1389https://doi.org/10.1126/science.aax8055
- 118.Subpopulations of vomeronasal sensory neurons with coordinated coexpression of type 2 vomeronasal receptor genes are differentially dependent on Vmn2r1European Journal of Neuroscience 47:887–900
- 119.Coordinated coexpression of two vomeronasal receptor V2R genes per neuron in the mouseMol Cell Neurosci 46:397–408https://doi.org/10.1016/j.mcn.2010.11.002
- 120.Evolution of spatially coexpressed families of type-2 vomeronasal receptors in rodentsGenome biology and evolution 7:272–285https://doi.org/10.1093/gbe/evu283
- 121.Loss of Kirrel family members alters glomerular structure and synapse numbers in the accessory olfactory bulbBrain Structure and Function 223:307–319
- 122.Kirrel3 is required for the coalescence of vomeronasal sensory neuron axons into glomeruli and for male-male aggressionDevelopment 140:2398–2408
- 123.A role for the EphA family in the topographic targeting of vomeronasal axonsDevelopment 128:895–906
- 124.Aberrant sensory innervation of the olfactory bulb in neuropilin-2 mutant miceJournal of Neuroscience 22:4025–4035
- 125.Neuropilin-2 mediates axonal fasciculation, zonal segregation, but not axonal convergence, of primary accessory olfactory neuronsNeuron 33:877–892
- 126.On the topographic targeting of basal vomeronasal axons through Slit-mediated chemorepulsionDevelopment 130:5073–5082
- 127.Robo-2 controls the segregation of a portion of basal vomeronasal sensory neuron axons to the posterior region of the accessory olfactory bulbJournal of Neuroscience 29:14211–14222
- 128.Olfactory sensory neuron-specific and sexually dimorphic expression of protocadherin 20Journal of Comparative Neurology 507:1076–1086
- 129.A role for TENM1 mutations in congenital general anosmiaClinical genetics 90:211–219
- 130.Differential requirements for semaphorin 3F and Slit-1 in axonal targeting, fasciculation, and segregation of olfactory sensory neuron projectionsJournal of Neuroscience 24:9087–9096
- 131.Molecular and structural basis of olfactory sensory neuron axon coalescence by Kirrel receptorsCell reports 37
- 132.Chemoreception regulates chemical access to mouse vomeronasal organ: role of solitary chemosensory cellsPLoS One 5https://doi.org/10.1371/journal.pone.0011924
- 133.Vomeronasal pump: significance for male hamster sexual behaviorScience 207:1224–1226
- 134.Access of large and nonvolatile molecules to the vomeronasal organ of mammals during social and feeding behaviorsJournal of chemical ecology 11:1147–1159https://doi.org/10.1007/BF01024105
- 135.Possible pheromone-carrier function of two lipocalin proteins in the vomeronasal organEMBO J 13:5835–5842https://doi.org/10.1002/j.1460-2075.1994.tb06927.x
- 136.How is the olfactory map formed and interpreted in the mammalian brain?Annual review of neuroscience 34:467–499https://doi.org/10.1146/annurev-neuro-112210-112917
- 137.Odorant receptor-derived cAMP signals direct axonal targetingScience 314:657–661https://doi.org/10.1126/science.1131794
- 138.Agonist-independent GPCR activity regulates anterior-posterior targeting of olfactory sensory neuronsCell 154:1314–1325https://doi.org/10.1016/j.cell.2013.08.033
- 139.A neuronal identity code for the odorant receptor-specific and activity-dependent axon sortingCell 127:1057–1069https://doi.org/10.1016/j.cell.2006.10.031
- 140.Imaging Neuronal Responses in Slice Preparations of Vomeronasal Organ Expressing a Genetically Encoded Calcium SensorJournal of Visualized Experiments
- 141.Near-optimal probabilistic RNA-seq quantificationNat Biotechnol 34:525–527https://doi.org/10.1038/nbt.3519
- 142.Modular, efficient and constant-memory single-cell RNA-seq preprocessingNat Biotechnol 39:813–818https://doi.org/10.1038/s41587-021-00870-2
- 143.Ensembl 2021Nucleic Acids Res 49:D884–D891https://doi.org/10.1093/nar/gkaa942
- 144.R: A language and environment for statistical computingR Foundation for Statistical Computing. Vienna, Austria
- 145.EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing dataGenome Biol 20https://doi.org/10.1186/s13059-019-1662-y
- 146.Integrated analysis of multimodal single-cell dataCell 184:3573–3587https://doi.org/10.1016/j.cell.2021.04.048
- 147.Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic DataCell Syst 8:281–291https://doi.org/10.1016/j.cels.2018.11.005
- 148.reticulate: Interface to ‘Python’
- 149.Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regressionGenome Biol 20https://doi.org/10.1186/s13059-019-1874-1
- 150.glmGamPoi: fitting Gamma-Poisson generalized linear models on single cell count dataBioinformatics 36:5701–5702https://doi.org/10.1093/bioinformatics/btaa1009
- 151.Clustering trees: a visualization for evaluating clusterings at multiple resolutionsGigascience 7
- 152.msigdbr: MSigDB gene sets for multiple organisms in a tidy data formatR package version 7
- 153.Fast gene set enrichment analysisbioRxiv 60012https://doi.org/10.1101/060012
- 154.Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profilesProc Natl Acad Sci U S A 102:15545–15550https://doi.org/10.1073/pnas.0506580102
- 155.Trajectory-based differential expression analysis for single-cell sequencing dataNat Commun 11https://doi.org/10.1038/s41467-020-14766-3
- 156.vegan: Community Ecology Package
- 157.GeneOverlap: Test and visualize gene overlaps
- 158.“Circlize” implements and enhances circular visualization in R
- 159.QuPath: Open source software for digital pathology image analysisSci Rep 7https://doi.org/10.1038/s41598-017-17204-5
Article and author information
Author information
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
Copyright
© 2024, Hills et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
- views
- 351
- downloads
- 13
- citations
- 0
Views, downloads and citations are aggregated across all versions of this paper published by eLife.