Introduction

The Vomeronasal organ (VNO), part of the vertebrate accessory olfactory system, is a major pheromone sensing organ, that is thought to be specialized to evoke innate social behaviors in mammals. The rodent VNO neuroepithelium consists of two major neuronal subtypes, that are defined primarily based on the expression of two distinct families of G-protein coupled receptors (GPCRs) and their associated G protein alpha subunit: vomeronasal type-I GPCRs (V1R) with G⍺i2 subunit (Gnai2) neurons or vomeronasal type-2 GPCRs (V2R) with G⍺o subunit (Gnao1) neurons. V1R and V2R expressing neurons are spatially segregated within the neuroepithelium, project to distinct locations in the accessory olfactory bulb along the anterior-posterior axis. Functional evidence suggests that V1R and V2R neurons may recognize ligands that mediate different behaviors (Chamero et al., 2011; Oboti et al., 2014; Trouillet et al., 2019; Pallé et al., 2020). However, we currently lack a comprehensive mapping of receptor-ligand interactions and how these lead to distinct behaviors. Cognate ligand binding in both Gnao1 and Gnai2 neurons, is thought to signal via the phospholipase-c pathway which results in the activation of the VNO specific cation channel, Trpc2 (Liman, Corey and Dulac, 1999; Stowers et al., 2002; Zufall, 2005)

In the main olfactory epithelium (MOE), neurons regenerate throughout life, continuously differentiating from a stem cell population into neurons and supporting cell types (Fletcher et al., 2017). Neuronal differentiation is also associated with the expression of a single GPCR to the exclusion of many others (Hanchate et al., 2015; Tan, Li and Xie, 2015), leading to the hallmark feature termed as “one neuron one receptor” (ONOR) rule (Malnic et al., 1999; Serizawa et al., 2000, 2003). VNO sensory neurons (VSNs) also continuously regenerate (Barber and Raisman, 1978; Wilson and Raisman, 1980) from a stem cell population, and differentiating cells take on the identity of supporting cells and mature neurons of either the Gnao1 or Gnai2 sub-types (Katreddi et al., 2022). While Gnai2 neurons seem to largely follow the ONOR rule in expressing one V1R gene per neuron (Serizawa, Miyamichi and Sakano, 2004), Gnao1 neurons seemingly deviate from this rule by co-expressing a combination of GPCRs. Thus, each Gnao1 neuron expresses a single V2R gene from the ABD family as well as a member of the C family V2R genes (Martini et al., 2001; Silvotti et al., 2007, 2011; Ishii and Mombaerts, 2011). In addition, further molecular diversification of Gnao1 neurons occurs due to combinatorial co-expression of a family of non-classical Major Histocompatibility Complex (MHC) genes, H2-Mv, that are co-expressed along with V2Rs (Ishii, Hirota and Mombaerts, 2003; Loconto et al., 2003; Ishii and Mombaerts, 2008). A few combinations of V2R and H2-Mv expression have been identified, indicating that these co-expression patterns are non-random and could be functionally important (Martini et al., 2001; Ishii, Hirota and Mombaerts, 2003; Loconto et al., 2003; Silvotti et al., 2007, 2011; Ishii and Mombaerts, 2008, 2011). However, identifying the precise combinatorial expressions of V2R-ABD, V2R-C and H2-Mv has been elusive, partly due to the challenges in generating RNA or antibody reagents to identify closely homologous genes, combined with the challenges with multiplexing more than three targets.

Starting with a common progenitor population, how do cells of the VNO sensory epithelium diversify in their gene expression patterns that define the sensory neuron types? Recent studies using single cell transcriptomic analysis have implicated Notch signaling along with downstream transcription factors, Tfap2e, Bcl11b to play a role in specifying the fates of developing VSNs (Enomoto et al., 2011; Katreddi et al., 2022; Lin et al., 2022). As the progenitor cells differentiate to mature neurons, transcriptional reprogramming downstream to Notch signaling leads to distinct neuronal populations: Gnao1 and Gnai2 neurons. While these studies identified transcription factors, it is not clear whether the two major differentiated VSN subtypes expressing Gnao1-V2Rs and Gnai2-V1Rs fundamentally differ in their cellular properties.

Here we applied a multilevel transcriptomics approach which bundles single cell RNA sequencing (scRNA seq) validated by low-level spatial transcriptomics to characterize gene expression of the mouse VNO sensory epithelium. Our results identified genes specifically expressed in sensory, non-sensory cell types and the divergence of gene expression between Gnao1 - Gnai2 neurons at different maturation stages. Analysis of mature neurons provided a comprehensive picture of specific co-expression patterns of VR and H2-Mv genes. Surprisingly, mature Gnao1 neurons revealed a significant enrichment of endoplasmic reticulum related gene transcripts and proteins, which is indicative of a fundamental divergence in the cellular function between these cell types. Electron microscopy further revealed a hypertrophic ER ultrastructure consisting of highly sinusoidal cubic membrane morphology, more prevalent in Gnao1 neurons.

Results

Single cell transcriptome of vomeronasal neuroepithelium

We standardized VNO neuroepithelium dissociation to obtain a single cell suspension along with optimizing cell viability. We performed single cell droplet-based capture coupled with scRNA seq, using the 10X genomics platform on four biological replicates consisting of male and female animals. Initial quality control steps were performed (see methods) to remove cells with high or low total RNA content and red blood cell contamination, resulting in the retention of 9185 cells from all samples for further analysis. Since our comparison of gene expression obtained from male versus female mice, did not reveal appreciable differences other than known sexually dimorphic genes that include Eif2s3y, Ddx3y, Uty and Kdm5d (Figure 1-figure supplement 1, supplementary table 2), we pooled cells from both genders for further analysis.

To group cells based on their gene expression profiles, we performed dimensionality reduction and unsupervised clustering, which resulted in 18 distinct clusters (Figure 1A), with gene expression markers specific to each cluster (Figure 1B). Major markers for each cluster were identified by differential expression analysis and the cell type identity was assigned to each cluster based on known markers for each cell type. Amongst 9185 cells, we observed that 63.7 percent (5852 cells) were neuronal cell types at different developmental stages which includes progenitor cells (Cluster 5; defined by expression of Ascl1, Neurod1, Neurog1), immature neurons (Cluster12; expression of Gap43), Gnao1 neurons (Cluster 14; based on expression of Gnao1, V2Rs), and Gnai2 neurons (Cluster 7, 16; characterized by Gnai2, V1Rs expression) (Figure 1C). The remaining 36.3 percent cells were classified to be of non-neuronal types which include sustentacular cells (Cluster 15; defined by Cyp2a5, Ppic), cells from non-sensory epithelium (Cluster 11; defined by Ly6d, Krt5), immune cells which include macrophages (Cluster 8, 6, 10: defined by C1qa, C1qb, H2-Aa), T cells (Cluster 9; defined by Cd3d, Cd3g), neutrophils (Cluster 2; defined by s100a8, s100a9), olfactory ensheathing cells (Cluster 3; defined by Gfap, Plp1, Sox10), endothelial cells (cluster 1 and 17; defined by Gng11, Egfl7, Eng), fibroblasts (Cluster 4; defined by Igfbp4, Dcn, Lum) and solitary chemosensory cells (Cluster 18; defined by Trpm5, Gnat3) along with a cluster that we could not attribute to a specific cell type (cluster 13; Obp2a, Gnai2 and Lcn3). A complete list of top 20 differentially expressed genes associated with each cluster and proportion of cells in each cluster are provided in supplementary table 1 and 3.

Figure 1.

We performed RNA in situ hybridization (RNA-ISH) on VNO sections to confirm the spatial classification of the various non-sensory cell type clusters based on some of the genes identified in our study (Figure 2A). Thus, we confirmed Ppic to be expressed in sustentacular cells and Ly6d in cells of non-sensory epithelium bordering the VNO lumen. Since, Cbr2 - a known marker of sustentacular cells is also expressed in non-sensory epithelium (Figure 1C), Ppic and Ly6d are better markers to distinguish them. Plp1 and Mgp, markers of olfactory ensheathing glial cells and endothelial cells/fibroblasts, respectively, were observed to be expressed in these cell types located at the periphery of the neuro-epithelium. Macrophage cell markers C1qb and H2-Aa were observed to be specifically expressed in scattered cells throughout the neuroepithelium, with H2-Aa co-expression also seen in sustentacular cells (Figure 2A), confirming the presence of these cell types within the neuroepithelium.

Figure 2.

Since our study identified multiple immune cell clusters, we examined these in some detail. Three clusters (Figure 1A, clusters 6, 8, 10 termed as Mac1, Mac2, Mac3 respectively) comprised of macrophage like cells based on expression of C1qa, C1qb, C1qc and MHC class II gene -H2-Aa (Figure 2B, 2C). Further analysis indicates that cluster 10 (Mac3) expresses classical microglial markers such as Tmem119, Aif1 (Iba1), P2ry12, Cx3cr1, Trem2, along with chemokine genes - Ccl3, Ccl4 and Ccl7 and Interleukins – Il1a, Il1b (Figure 2C), indicating an activated macrophage subset. Cluster 8 (Mac2) is similar to Mac3 in regard to expression of classical tissue resident macrophage marker - Adgre1 (F4/80), C1qa, C1qb, C1qc but lacks the expression of chemokine and Interleukin genes mentioned earlier. In comparison to Mac3 and Mac2, the Mac 1 cluster is distinguished by the lack of C1qa, C1qb, Cd68, Adgre1(F4/80), microglial markers but possesses high levels of Napsa, Lsp1, H2-Aa (Figure 2C, Supplementary table 1, 4). In summary, the identification of discrete macrophage cell types from our data, along with experimental validation of their presence using some of the expression markers, points to the possibility of tissue resident or tissue-infiltrating macrophages within the VNO neuroepithelium. These VNO macrophage types could have functional implications towards influencing processes such as neuronal differentiation or responses to tissue repair and injury.

On the other hand, we did not observe evidence via RNA-ISH for the T-cell markers (Cd3d, Cd3g), raising the possibility that these immune cell types may have co-purified from blood. Lastly, we were able to identify via RNA-ISH, the expression of Trpm5 (Figure 2A) and Gnat3 genes from our dataset, that mark solitary chemosensory cells expressing taste receptors and related markers (Ogura et al., 2010), indicating that rare cell-populations were captured, and their markers identified in our study.

Transient gene expression in developing neurons

Our next step was to perform in depth analysis of the neuronal cell types. We pooled cells from clusters representing neuronal types of different developmental stages (clusters 5, 7, 12, 14, 16) and created a new Seurat object and performed re-clustering. This resulted in 13 clusters (labelled n1-n13), which potentially represent distinct neuronal sub-populations within the broadly defined Gnao1 and Gnai2 neurons (Figure 3A, Figure 3-figure supplement 1A). Since neurogenesis is a continuous process in VNO neurons (Barber and Raisman, 1978; Wilson and Raisman, 1980), these clusters also represent a snapshot of cells at different developmental states, identified by known markers: progenitor cells (Cluster n5; defined by expression of Neurod1, Neurog1), immature neurons (Cluster n6, n7; defined by Gap43) expressing Gnao1, Gnai2. Among the mature neuronal clusters, four clusters (n1-n4) correspond to Gnao1 neurons, and 6 clusters (n8-n13) correspond to Gnai2 neurons (Figure 3–figure supplement 1B). Since immature neurons separated into two clusters expressing Gnao1 or Gnai2, these imply a possible branch point in the developmental trajectory. To confirm this, we performed pseudo time analysis on this object using Slingshot, and the progenitor neurons (cluster n5) as a start/anchor point. The trajectory of development confirms a split of immature Gap43+ neurons into Gnao1 or Gnai2 cell clusters (Figure 3A). We performed differential expression analysis between Gap43+; Gnao1+ (cluster n6) and Gap43+; Gnai2+ (cluster n7) immature neurons (Figure 3B, Supplementary Table 5). To identify potential transcription factors in this developmental branch point, we manually curated genes that have molecular features associated with DNA binding domains or transcription regulation based on literature (Figure 3C).

Figure 3.

To our surprise, we found multiple genes that are known transcription regulators but are not previously reported to be enriched in Gnao1 neurons. These include Creb5, Prrxl1, Shisa8, Lmo4 and Foxo1. Creb5, Prrxl1 are transcription factors expressed only in immature Gnao1 neurons whereas Shisa8, Lmo4 and Foxo1 expression is enriched in immature Gnao1 neurons compared to mature ones (Figure 3B). Prrxl1 was reported to auto-repress its expression which explains how its expression is limited to immature Gnao1 neurons (Monteiro et al., 2014).These data indicate that the developmental transition of VNO neurons to Gnao1 lineage may involve the transient expression of transcription factors within immature neurons after dichotomy is established.

Divergent gene expression profile of mature Gnao1 and Gnai2 neurons

Having identified differences in gene expression of immature Gap43+; Gnai2+ and Gap43+; Gnao1+ neurons, we focused on mature Gnao1 and Gnai2 neurons. We performed differential gene expression analysis of mature Gnao1 (cluster 14 from Figure 1) and Gnai2 neurons (cluster 7, 16 from Figure 1) and identified 924 differentially expressed genes with at least log2 fold change of 1 and p-value<10e-6. Of these, 456 genes are enriched in Gnao1 neurons and 468 genes in Gnai2 neurons (Figure 4A). In agreement with previous studies, we see enrichment of V2Rs, H2-Mvs, B2m (Ishii, Hirota and Mombaerts, 2003; Loconto et al., 2003), Robo2 (Prince et al., 2009), Tfap2e (Lin, Jennifer M et al., 2018), Calr4 (Dey and Matsunami, 2011) in Gnao1 neurons and V1Rs, Meis2 (Chang & Parrilla, 2016) in Gnai2 neurons (Figure 4A, Supplementary Table 6). Along with such known genes, our analysis revealed new differentially expressed genes in the two neuronal subtypes.

Figure 4.

To confirm the differential gene expression amongst mature neurons, we performed two-color fluorescence RNA-ISH of the top candidates from Figure 4A, using Gnao1 or Gnai2 as co-labelling markers (Figure 4B, 4C, 4D). Among these genes, Nrp2 and Socs2 expression is restricted to Gnai2 neurons (Figure 4C); Tppp3 and Gm36028 expression is restricted to Gnao1 neurons (Figure 4D). Some of the enriched genes such as Ckb in Gnai2 and Krt18 in Gnao1 neurons, were also observed to be expressed in sustentacular cells (Figure 4C, 4D). Furthermore, from chromogenic RNA-ISH experiments, we observed genes such as Nsg1, Rtp1, Dner, Qpct, Pcdh7 to be either restricted to or highly expressed in Gnai2 neurons (Figure 4–figure supplement 1A) and others such as Apmap, Selenof, Hspa5, Itm2b, Agpat5 were observed to have a gradient of high to low level expression from Gnao1 to Gnai2 zones (Figure 4–figure supplement 1B). Two genes: Sncg and Prph specifically expressed in scattered subsets of Gnao1 or Gnai2 neurons respectively (Figure 4–figure supplement 1A, 1B) indicating selective expression in a particular neuronal subpopulation.

These experiments validated that several of the candidate genes identified as differentially expressed from the scRNA seq data, are highly specific and spatially localized to one of the neuronal sub-populations. A complete list of differentially expressed genes is in supplementary table 6. Taken together our data indicate that VNO neurons that start from a common progenitor, go on to transiently express genes during immature stages that could decide their cell fate and these cells exhibit markedly different gene expression profiles in the mature stage.

Differential ER environment in Gnao1 and Gnai2 neurons

Given that Gnao1 and Gnai2 neurons express divergent GPCR families, we asked the question whether the differentially expressed genes indicate fundamental differences in cellular or molecular processes between these two cell types. We performed gene ontology (GO) enrichment analysis on the two sets of genes. Interestingly, we observed that amongst the Gnao1 neuron genes, most of the enriched terms with p value <0.05, are related to endoplasmic reticulum (ER) processes, which include ‘protein folding’, ‘response to ER stress’ and ‘ERAD pathway’ (Figure 4E). Conversely, a gene ontology analysis of Gnai2 neuron genes did not reveal any particular enrichment of GO terms in this neuronal subset. As more than 20% of the genes associated with protein folding (Figure 4E) are enriched in Gnao1 neurons, we sought to probe the spatial transcription pattern of some of the candidates. Fluorescence RNA in-situ hybridization confirmed the selective expression of ER genes such as Creld2, Dnajc3, Pdia6 and Sdf2l1 amongst Gnao1 neurons (Figure 5A). To test whether these findings extend to the protein level, we performed immune labelling of VNO sections using antibodies specific to a collection of ER proteins (Figure 5B, 5C). Surprisingly, several of the antibodies such as anti-Hspa5, anti-PDI, anti-Atlastin1, anti-NogoB and anti-SEKDEL detected ER protein localization selectively amongst the basal zone Gnao1 neurons, while some others such as anti-Sec61β, anti-Calnexin, anti-Grp94, anti-Reep5, anti-Climp63 detected a pattern of Gnao1 neuron enrichment, along with lower levels in Gnai2 neurons. Amongst the ER proteins tested by immune labelling, we observed that RNA level of Atlastin1 and PDI were comparable across Gnao1, Gnai2 neurons in our dataset (Figure 5–figure supplement 1), indicating that enrichment of ER proteins in Gnao1 neurons could also arise due to post-transcriptional regulatory mechanisms.

Figure 5.

Many of the Gnao1 enriched ER genes and their protein products are chaperones, while others, such as Atlastin1, Reep5, NogoB and Ckap4/Climp63 are ER structural proteins (Goyal and Blackstone, 2013). The upregulation of both types of ER proteins prompted us to examine the cells by electron microscopy. Serial sections (70 nm thick) of ∼1 mm square regions of the VNO were collected on tape and were imaged by scanning electron microscopy (Baena et al., 2019). Since the sections contained the entire VNO, we were able to distinguish Gnai2 from Gnao1 cells by their location as they are spatially segregated into apical or basal zones respectively within the neuroepithelium (Figure 6A). In Gnao1 cells, more than half of the cytoplasm contained dense smooth ER with cubic membrane morphology (Almsherqi, Kohlwein and Deng, 2006; Almsherqi et al., 2009), also known as organized smooth ER (OSER) (Snapp et al., 2003) (Figure 6A, 6B, Figure 6–figure supplement 1). In the remainder of the cytoplasm, there were closely apposed parallel sheets of ER membrane that were contiguous with lamellar ER surrounding the nucleus (Figure 6–figure supplement 1). The apical/Gnai2 neuronal cell bodies are smaller than those of Gnao1 cell bodies. They also appeared to have cubic ER membranes, but these seemed denser and occupied a smaller fraction of the cytoplasmic volume (Figure 6A, 6C).

Figure 6.

Taken together, our data indicates that mature Gnao1 neurons differ substantially from Gnai2 neurons in selectively upregulating the transcription of several ER genes. Furthermore, this neuronal subset also exhibits higher levels of ER proteins, potentially due to post-transcriptional mechanisms, as well as a hypertrophic smooth ER that is arranged in the form of organized cubic membrane ultrastructure. Collectively, these observations indicate a specialized ER environment that could play a significant functional role in the Gnao1 subset of VNO neurons.

Co-expression of Vomeronasal receptors

One of the differentiating features of vomeronasal sensory neurons from the main olfactory sensory neurons is the deviation from ‘one neuron one receptor rule’, especially in mature Gnao1 neurons. V2R receptors have been classified into families-A, B, D, E and a distinct family-C that is phylogenetically closer to calcium sensing receptor (Young and Trask, 2007). Several studies have established a broad pattern of co-expression amongst these V2R family members, such as members of the family ABDE co-express with family-C GPCRs and components of the family-C (Vmn2r1 -also termed -C1, Vmn2r2-Vmn2r7 -termed C2) in turn co-express with each other in specific combinations (Martini et al., 2001; Silvotti et al., 2007). Furthermore, H2-Mv genes within Gnao1 neurons also exhibit non-random co-expression patterns with the V2R GPCRs. (Ishii, Hirota and Mombaerts, 2003; Loconto et al., 2003; Silvotti et al., 2007, 2011; Ishii and Mombaerts, 2008, 2011). To further investigate these co-expression patterns and subpopulations within Gnao1 neurons, we took a closer look at this neuronal cell type, within our scRNA seq data. From the unsupervised clustering of neuronal cell types obtained from Figure 3A, Figure 3 supplement A, we separated mature Gnao1 neurons (cluster n1-n4) for further analysis. These four clusters (Figure 7A) are found to be organized majorly based on the expression of family-C V2Rs and H2-Mv genes, especially Vmn2r1 or Vmn2r2 (Figure 7B-E). Of the 4 clusters, cluster n1 and cluster n4 express Vmn2r1 (Figure 7B, 6E) whereas cluster n3 neurons are distinguished by the co-expression of Vmn2r2 along with other family-C members (Vmn2r3, Vmn2r4, Vmn2r5, Vmn2r6. Vmn2r7) (Figure 7C, 6E). Cluster n2 is a mixture of either Vmn2r1 or Vmn2r2 expressing neurons (Figure 7E), with Vmn2r2 being similarly co-expressed along with other family-C members in those cells. We observed that among family-ABDE members, A1-A6 V2Rs are highly expressed in cluster n2, n3 which co-express Vmn2r2 whereas A8-A10, B, D, E family V2Rs are mainly expressed along with Vmn2r1 in cluster n1 neurons (Figure 7E). These observations of family biased expression of ABDE with Vmn2r1 or Vmn2r2 agrees with the earlier studies (Silvotti et al., 2007, 2011) that have used antibody or RNA probe based approaches.

Figure 7.

Furthermore, we observed that the expression of H2-Mv genes (H2-M10.1, M10.2, M10.3, M10.4, M10.5, M10.6) is largely confined to clusters n2 and n3, where they are co-expressed selectively in Vmn2r2 positive neurons and mostly absent from Vmn2r1 expressing cluster n1, n2 neurons (Figure 7D, 7E, Figure 7–figure supplement 2A). These observations align with prior reports from antibody or RNA probe-based experiments, that demonstrated a sub-division amongst Gnao1 neurons based on their expression of family-C Vmn2rs and the co-expression of H2-Mv genes with the Vmn2r2 subset (Silvotti et al., 2007).

We then asked whether at the single neuron level, we could identify trends in co-expression patterns between ABD-family V2Rs and H2-Mv genes, or between members of the ABD family themselves. We performed co-expression analysis to identify specific VR and H2-Mv combinations by setting an expression threshold from normalized counts. The distribution of normalized expression for VRs and H2-Mvs across cells was bimodal, with the first peak near zero and a second peak likely corresponding to true expression (Figure 7–figure supplement 1D-F). We set the starting of the second peak as a threshold to call the expression of VRs or H2-Mvs in a single cell. Thus, a VR or H2-Mv was considered as co-expressed in a cell only if the normalized expression value surpasses the set threshold based on the distribution (Figure 7–figure supplement 1E, 1F). The threshold was 1.25 for V2Rs (Figure 7–figure supplement 1E) and H2-Mv genes (Figure 7–figure supplement 1F). This analysis resulted in identification of multiple co-expressing V2Rs and H2-Mv genes per cell. The chance that a combinatorial pattern identified can be due to a doublet is dependent on the frequency of expression of that gene. Therefore, for abundantly expressed family C V2Rs, it is possible some combinatorial patterns may be from doublets. But in the case of family ABD V2Rs, the probability of a particular combination falsely identified due to a doublet more than once is very low. Hence, we assigned a threshold for combinatorial co-expression patterns that are identified in 2 or more cells to be more reliable. These combinatorial expression patterns along with their observed frequency are listed in supplementary table 7.

Within the expression patterns represented by a frequency of more than two cells, we observed that the family C1 V2R (Vmn2r1) is mostly expressed alone without any other family-C members (Figure 7F, 7E) as shown earlier (Silvotti et al., 2007). Surprisingly a few cells (31 neurons out of 1025 Gnao1 neurons) were observed to co-express Vmn2r1 along with Vmn2r7 (Supplementary table 7), contrary to earlier reports that Vmn2r1 does not co-express with family-C2 receptors. Most of Gnao1 neurons with family C2 V2Rs (Vmn2r2-2r7) co-express multiple members amongst them, ranging from 2 to 6 V2Rs per cell (Figure 7F, 7E).

In the case of family-ABD V2R expression, we observed that most of the Gnao1 neurons follow ‘one ABD-V2R, one neuron’ rule, except for few that express two ABD-V2Rs (Figure 7G). Amongst cells that express more than two ABD-V2Rs, the combination of Vmn2r20 and Vmn2r22 stood out as the highest (44 cells), followed by Vmn2r39, Vmn2r43 and Vmn2r50 (10 cells, Supplementary table 7). It is worth noting that these ABD co-expressions in each cell are observed irrespective of whether that cell expresses famly-C1 or C2 V2R (Figure 7G). Among the Gnao1 neurons (n1-n4), 229 cells did not express family-C V2Rs (Figure 7F) and 46 cells did not express ABD-V2Rs (Figure 7G) as per our threshold.

In the case of H2-Mv genes, each cell expresses 1-4 members mostly with famlily-C2 V2Rs (Figure 7H, 7D). Interestingly, in deviation from this, we observed a pattern whereby H2-M1, -M9 and -M11 genes that show sequence divergence from the rest of the H2-Mv cluster (Ishii, Hirota and Mombaerts, 2003), are specifically expressed in cluster n4 and restricted cells amongst cluster n3, where they co-express with Vmn2r1 GPCR, rather than family-C2 V2Rs (Figure 7–figure supplement 2B, 2D). Overall, 56.9% of Gnao1 cells that express either C1-V2R or C2-V2R, did not express any H2-Mv genes (Figure 7H) as per our threshold, reinforcing earlier observations (Ishii and Mombaerts, 2008) that their functional role maybe restricted to a subset of V2R expressing neurons. Furthermore, these H2-M1, H2-M9 or H2-M11 expressing cells along with Vmn2r1, also co-express either Vmn2r81 or Vmn2r82 (V2rf) (Figure 7–figure supplement 2B, 2C) as identified before (Ishii, Hirota and Mombaerts, 2003). Two color RNA in situ hybridization with a common probe targeting Vmn2r81 and 82 and separate probes for H2-M1, H2-M9 and H2-M11 shows that Vmn2r81/82 cells always colocalized with the selected H2-Mv probes (Figure 7–figure supplement 2E). Overall, these data validate the co-expression analysis methodology and identify novel combinatorial co-expression patterns of V2Rs and H2Mv genes.

We next asked the question whether the total GPCR or H2-Mv expression level is regulated or capped in a single cell? In other words, when a cell expresses multiple V2Rs, does their combined expression remain similar to cells that express one V2R or H2-Mv gene? To answer this question, we calculated the total family-C-V2R, ABD-V2R and H2-Mv expression in a single cell that either express one member or multiple members of each family (Figure 7I-7K). The results showed that when cells express more than one V2R, the total V2R expression of either family-C or family-ABD (Figure 7I-7J) scales up proportional to total number of expressed V2Rs of that family. This is also true in the case of multiple H2-Mv expression where, total H2-Mv expression is proportional to the number of H2-Mv genes expression in a cell (Figure 7K). As a control, we see that the levels of Gnao1 does not change with number of V2Rs or H2-Mv genes expressed in a cell (Figure 7I-7K) indicating that these measurements are form singlet cells.

Lastly, we performed a similar co-expression analysis for V1Rs with a threshold value of 2.5 (Figure 7–figure supplement 1D). In similarity to the pattern observed for ABD-V2Rs, most Gnai2 neurons expressed a single V1R per cell with some cells co-expressing two V1Rs (Figure 7–figure supplement 3B, 3A, Supplementary table 8). Some of these V1R combinations have been reported recently (Lee, Kume and Holy, 2019). Even in the case of Gnai2 neurons, total V1R expression in the cells that co-express two receptors is higher compared to cells that express one GPCR with Gnai2 levels being same (Figure 7–figure supplement 3C). Collectively, these data reveal the patterns of GPCR, H2-Mv co-expression in VNO neurons, which is likely to be instrumental in deciphering how these neurons respond to single or combination of stimuli, as well as the cellular mechanisms that orchestrate these expression patterns.

Discussion

The vomeronasal organ has been a model of considerable interest for molecular sensory biology, to study the diversification of sensory neurons, uniquely evolved gene families and their patterns of expression, all of which are essential to understand how specific social chemo-signals elicit innate behaviors. In this study, using single cell and low-level spatial transcriptomics, we developed a comprehensive resource identifying cell types of the VNO neuroepithelium and studied the developmental and receptor co-expression patterns within sensory neurons. We organized our dataset as an interactive web portal resource, accessible from www.scvnoexplorer.com, where a gene can be queried for its expression pattern and levels amongst cell clusters.

The identification of distinct cell clusters and the validation of marker expression through RNA in situ hybridization, revealed a heterogeneous cellular composition of the vomeronasal neuroepithelium. In particular, we identified novel and specific gene markers, such as Ppic in sustentacular cells and Ly6d in non-sensory epithelium, which would be helpful in future to genetically tag or target these cell types and to study their developmental origins. Our study also identified specific types of macrophage cells within the neuroepithelium. Given the increasing appreciation of the role of tissue resident or infiltrating macrophages in a variety of physiological processes (Nobs and Kopf, 2021; Mass et al., 2023) it is tempting to speculate that macrophages could also play an important role in VNO neuronal differentiation or the regulation of inflammatory/immune responses to external insults. Future experiments, performing targeted purification, deep sequencing, and functional manipulation of VNO macrophage types could provide deeper insights on their role in VNO function and whether they share cellular identity with MOE macrophages.

Development of VSNs

VNO sensory neurons expressing Gnao1 and Gnai2 develop continuously throughout life from a single progenitor cell population located within the VNO marginal zones and differentiate to immature and mature neurons. Recent developmental studies have demonstrated that the dichotomy is established by Notch signaling followed by the transcription factor, Bcl11b expression that marks the progenitor cells to take the Gnao1 fate (Enomoto et al., 2011; Katreddi et al., 2022). Most Importantly, a critical period of Notch signaling determines whether the progenitors give rise to apical, basal or sustentacular cell types. After the Gnao1/Gnai2 identity is established, the transcription factor, Tfap2e, is continuously expressed starting from Gap43+, Gnao1+ neurons and persists into mature Gnao1 neurons and has been shown to be required to maintain the Gnao1 fate (Lin et al., 2018, 2022; Katreddi et al., 2022). Pseudo-time analysis on our data, confirmed a developmental trajectory that starts with progenitor cells, leading to immature neurons that take on either Gnai2 or Gnao1 fate. From our scRNA seq data, a comparison of these immature neurons revealed that apart from Tfap2e, transcription factors: Creb5, Prrxl1, Lmo4, Foxo1 also express at the same time during the immature stage of Gnao1 neurons. But unlike Tfap2e, the expression of Creb5 and Prrxl1 was observed to be restricted to immature Gnao1 neurons. Although, Foxo1 and Lmo4 continue to be expressed in mature Gnao1 neurons, their RNA levels are high in immature Gnao1 neurons and reduce upon maturation. Among the transcription factors identified, Prrxl1 has been shown to involved in the development of dorsal root ganglion neurons where it has been shown to autoregulate its expression (Monteiro et al., 2014)

Our observations on the tight temporal regulation of these transcription factor related genes in the developmental trajectory of Gnao1 neurons indicates that additional molecular players might be important in further specification or maintenance of this neuronal lineage. Further experiments are needed to tease out the mechanism and precise choreography by which these transcription factors collectively specify and maintain the Gnao1 cell identity as well as the gene expression patterns unique to these neurons.

Receptor co-expression in Gnao1 neurons

One of the key features distinguishing basal zone Gnao1 VNO neurons from the apical Gnai2 or MOE neurons, is the combinatorial co-expression of two or more V2R family GPCRs along with the H2-Mv class-1b MHC molecules in each cell. Our analysis of receptor expression in mature VNO neurons, revealed that members of the V1R and V2R-ABD gene families largely follow the one-neuron-one receptor rule, but with multiple exceptions. Some V1Rs were found to consistently co-express in Gnai2 neurons and the major combinations observed from our analysis (such as Vmn1r85/86, Figure 7–figure supplement 3A), match those identified recently from a functional single cell study (Lee, Kume and Holy, 2019). leading confidence to the methodology used in our co-expression analysis. Receptors of the ABD-V2R subfamilies are also mostly expressed at one receptor per neuron, with notable exceptions being the observed combinations of Vmn2r20/22 and Vmn2r39/43/50 being co-expressed per cell.

In our study, the analysis of family-C, family ABD-V2Rs and H2-Mv genes expressed per cell, confirmed the earlier reported sub-division of the basal/Gnao1 neurons into those that express either Vmn2r1 (family-C1) or Vmn2r2-Vmn2r7 (family-C2), resulting in a broad four-part division of the mature VNO neurons: those expressing a) Gnai2/V1Rs, b) Gnao1/Vmn2r1, c) Gnao1/Vmn2r2-r7 with H2-Mv and d) Gnao1/Vmn2r2-r7 without H2-Mv. The co-expression of ABD-V2R sub-families A1-A5 and most of the H2-Mv genes with Vmn2r2/family-C2 neurons (Figure 7C-E, cluster n3, n2, Supplementary table-7) as well as the preferential co-expression of A6-A9, B, D, E sub-families with Vmn2r1 along with exclusion of H2Mv genes from these neurons (Figure 7E), points to the overall non-random logic of V2R/H2Mv co-expression. It would be interesting to see how receptor de-orphanization and future functional experiments will map onto these elaborate co-expression patterns. At the same time there were notable deviations observed from these rules, namely the restricted and sparse expression of H2-M1/M9/M11 with Vmn2r1 neurons in cluster n4 (Figure 7E, Figure 7–figure supplement 2B, 2E), where the receptors Vmn2r81/82 are also expressed selectively in these neurons.

Of note, the combined expression of VR/H2-Mv in a cell is not capped and is proportional to number of VRs and H2-Mv genes expressed in that cell (Figure 7I-7K). It would be interesting to see whether the levels of co-expressing transcripts also translate to protein levels in cells and their impact on functional interactions, if any, between co-expressed proteins.

ER environment in VNO neurons

When comparing gene expression amongst Gnai2/Gnao1 mature neurons, what stood out was the unexpected enrichment of gene ontology terms associated with ER functions within Gnao1 neurons (Figure 4A, E). Most of these ER-Gnao1 enriched genes from our dataset, match the ones identified from a similar comparison from a recent study by Lin et al 2022, supporting our findings. We validated several of these ER genes via RNA-ISH to be either exclusively expressed or enriched in Gnao1 neurons. Even more puzzling was our observation that this enrichment pattern extended at the protein level to many ER proteins probed by immune labelling (Figure 5). In some cases, both RNA and proteins were Gnao1 enriched, whereas in others, the RNA levels were comparable across Gnai2/Gnao1 neurons, but the protein levels were biased towards Gnao1 cells. The enrichment of ER genes/proteins via transcriptional and/or post-transcriptional mechanisms presents a new feature of these cell types that could be associated with their differentiation and function. High levels of several chaperones such as Creld2, Dnajc3, Pdia6, PDI, Grp78/Bip, Grp94, Calnexin in Gnao1 neurons (Figure 5B) compared to Gnai2 neurons indicates a chaperone rich ER environment in these neurons. The chaperone Calr4 has been shown to be enriched in Gnao1 neurons (Dey and Matsunami, 2011). Our observations indicate that the enrichment of chaperones is a generalized feature of Gnao1 neurons that extends to several other genes/proteins of related function and localization.

Since it has been proposed that the V2R family is an evolutionarily recent acquisition in rodents and marsupials (Takigami et al., 2004) and sub-family members, such as sub family-C2, A1-A6 V2Rs, as well as H2-Mv genes have evolved in Muridae (Young and Trask, 2007; Silvotti et al., 2011), it is possible that mouse Gnao1 neurons may have required to co-evolve ER protein folding mechanisms to handle multiple GPCR/H2-Mv co-expression and their putative protein interactions. In addition to our observations on the generalized upregulation of ER genes in Gnao1 neurons, it may be possible that some specific ER genes are also necessary for proper folding of vomeronasal type-2 GCPRs. These data and observations assume significance given the well-recognized fact that the functional reconstitution of these GPCRs and H2-Mv molecules in heterologous cells remains a persistent challenge, with some success being reported using the chaperone Calr4 (Dey and Matsunami, 2011). On the other hand, olfactory receptors have been shown to co-opt the ER unfolded protein response (UPR) pathway via the ER proteins PERK/CHOP and transcription factor ATF5 as a mechanism for their functional expression, setting the stage for ONOR expression pattern (Dalton, Lyons and Lomvardas, 2013). This has been further extended to show that differential levels of ER chaperones help in axonal guidance of olfactory sensory neurons (Shayya et al., 2022). While ATF5 is expressed broadly in the VNO, across Gnai2 and Gnao1 neurons, further experiments are required to evaluate whether V2Rs and Gnao1 neurons adopt similar mechanisms and how the chaperone rich ER environment observed here might impact the expression, folding of GPCRs and the function of these neurons.

It is worth noting that in addition to the expression of ER chaperones in Gnao1 neurons, we also observed distinctly higher levels of ER structural and membrane proteins. Levels of the membrane curvature stabilizing proteins - Reep5 and NogoB/Reticulon4B, the three-way junction formation protein - Atlastin1 (Goyal and Blackstone, 2013), ER sheet lumen spacer protein - Climp-63/Ckap4 as well as Sec61β - a commonly used ER membrane marker were observed to be high in Gnao1 neurons (Figure 5B, 5C). Thus, it is possible that an added layer of complexity could involve modulation of the total ER content or ER structure in Gnao1 neuronal subtypes. Our electron microcopy results (Figure 6, Figure 6-figure supplement 1) revealed a strikingly higher smooth ER membrane content in Gnao1 neurons compared to Gnai2 neurons. Interestingly, ER in both neuronal types adopt a cubic membrane architecture, that is more pronounced in Gnao1 neurons, packing the cell body with sinusoidal membrane. Early electron microscopy studies in rodents, identified that smooth ER content is higher in VSNs than olfactory sensory neurons (Ciges et al., 1977; Breipohl et al., 1981; Taniguchi and Mochizuki, 1982), however, these studies preceded the discovery of Gnao1-Gnai2 dichotomy. To our knowledge, Breipohl and colleagues (Breipohl et al., 1981) reported the presence of smooth ER whorls in VSNs, later termed as cubic membranes. Future studies using serial section EM would be instrumental in reconstructing the 3D architecture of VSN ER.

Cubic membranes, representing highly curved periodic structures, have been observed in a variety of biological contexts, however their functional significance has been harder to pin down (Almsherqi, Kohlwein and Deng, 2006; Almsherqi et al., 2009). In the context of VNO, and in particular Gnao1 neurons, we speculate that the high ER content and dense packing in the form of cubic membranes, could indicate a perpetual stress like state in these neurons, necessary to address the unique folding requirements of V2R or H2-Mv proteins. For instance, the maturation of these proteins and their proper folding or multimerization may require a slower transit through the ER, necessitating the sinusoidal arrangement of ER membranes. Likewise, the V2R biosynthetic processes in Gnao1 neurons may require enhanced levels of ER chaperones, which could in turn induce the ER membranes into a homeostatic response. Can V2Rs themselves trigger the upregulation of ER chaperones and the hypertrophic ER, or are these cellular characteristics established before the onset of V2R expression? From our single cell transcriptomics data, we tried to identify at what stage ER chaperones express during the developmental trajectory of VSNs. We observed that onset of Vmn2r1 expression coincides with the upregulation of ER chaperones in cluster n6 immature Gnao1 neurons (Figure 8 A, B). However, some Gnao1 cells show upregulated ER gene expression without Vmn2r1/Vmn2r2 expression (Figure 8B), hence further experiments are needed to dissect the mechanism and precise relationship between neuronal differentiation, GPCR expression and ER chaperones as well as ER ultrastructure.

Figure 8.

Overall, our data opens a new aspect to look at Gnao1 neurons, when understanding their function: a highly specialized ER environment that is divergent from Gnai2 neurons. The ER genes and their Gnao1 biased expression we identified, could be downstream targets of specific transcription factors and as a result of neuronal differentiation, while functionally these could be required for the proper folding and co-expression of the V2R GPCRs. The developmental trajectory of VSNs into Gnai2/Gnao1 neurons via transcription factors and their differential ER environment is summarized as a model in Figure 8C.

In conclusion, the comprehensive scRNA seq analysis presented in this study contributes valuable insights into the complexity of the vomeronasal neuroepithelium, offering a roadmap for further investigations into the molecular mechanisms underlying sensory perception and neural development in this specialized sensory organ. Vomeronasal neurons may also serve as model system to study specialized ER structure-function relationship and its role in maturation of GPCR expressing neurons.

Methods

Animals

C57BL/6J mice were purchased from JAX (Strain #000664) and were used for all experiments. Mice were housed in a specific pathogen-free barrier facility, with a 12-hour light-dark cycle and ad-libitum provision of feed and water. For single-cell RNA sequencing experiments, male and female mice were weaned at postnatal 3 weeks age, followed by housing in separate individually ventilated cage racks to avoid exposure to opposite sex stimuli, and used at 7-8 weeks age. All experiments were carried out with approval from the Institutional Animal Ethics Committee of TIFR Hyderabad.

VNO dissociation

Papain dissociation buffer (PDB) was made by reconstituting single use papain vial from Worthington Biochemical Corporation with 5 mL Earle’s balanced salt solution (EBSS) (pH 7.2) and warmed to 37°C until solution appears clear by maintaining 95% CO₂, 5% O₂ environment. VNOs were dissected from male and female animals (6 male, 10 female) and processed separately. The sensory epithelium was separated using forceps and immediately placed in EBSS equilibrated with 95% CO_2, 5% O₂. 60U of DNase was added to the prewarmed PDB made earlier and 3mL of it is transferred to a single well of 12-well dish. Neuroepithelial tissue from multiple animals was transferred to this well containing the final dissociation buffer and was cut into small pieces and the suspension was transferred to 14 mL tube and was incubated with gentle shaking at 37°C for 30 minutes on a thermomixer (Eppendorf) by passing 95% CO₂, 5% O₂ on top of the liquid headspace. During this incubation, the solution was triturated with fire-polished Pasteur pipette every 10 minutes and at the end the incubation. After incubation, the suspension was passed sequentially through 70µm, 35µm filters to remove tissue debris. The filtered suspension was further layered over a mix of ovomucoid inhibitor (2.5 mg/mL) and albumin (2.5 mg/mL) in EBSS and centrifuged at 400 x g for 5 minutes at room temperature to remove subcellular debris or membrane fragments. The supernatant was discarded, and the cell pellet was used for next steps after Hank’s balanced Salt solution (HBSS) wash.

Single-cell library preparation and sequencing

Dissociated neurons resuspended in HBSS were used for library preparation. Using hemocytometer and trypan blue staining, cell density was estimated at 1200 cells/μl and 1100 cell/μl for male and female samples, respectively with greater than 80% cell viability. Single-cell capture, library preparation was done using a Chromium Next GEM Single Cell 3’ GEM, Library and Gel Bead Kit v3.1 on a 10X Chromium controller (10X Genomics). The volume of suspension to be loaded was decided as per the manufacturer recommendation to ensure the target capture of 6000 cells per well. Single cell suspension from male and female samples was loaded onto two separate wells of different chips giving a total of four libraries (2 male and 2 female samples). Each library was barcoded and sequenced separately on a single lane of HiSeq X to obtain a mean depth of at least 100,000 reads per cell in 2 x150bp configuration.

scRNA seq data analysis

Reference details of all software used for data analysis are mentioned in a tabular form below. Raw reads for each sample were trimmed as per requirements specified by 10X genomics (read 1 - 28 bps and read 2 - 91 bps) using trim galore. Trimmed reads were aligned separately to reference genome mm10 (Mus_musculus.GRCm38.dna.primary_assembly.fa) and the count matrix was generated using cellranger-5.0.1 for each sample. To avoid alignment errors, read-through transcripts were removed from the reference. Male (2) and female (2) replicates were integrated by using the cellranger’s ‘aggr’ option. Before downstream analysis stringent filtering was implemented to retain high quality cells as per the following. Ambient RNA contamination was removed using soupX software with default parameters in auto mode. The adjusted count matrix from soupX consisting of 10615 cells was used as an input to Seurat for further analysis. To remove potential doublets and low gene count cells, an additional filter was applied to select cells that express 200-7000 genes resulting in dropping 683 cells. Principal component analysis (pca) was performed, and the number of dimensions required for clustering was determined as 37 based on elbow plot and Jack Straw plot. Clustering was performed using Seurat’s FindClusters function with resolution=0.3. Markers for each cluster were identified by FindAllMarkers function. Cluster identity was assigned based on known markers of each cell type. To clean up the dataset further, cells (408) expressing Hbb-bs gene were considered as RBC contamination and a cluster (215 cells) enriched with mitochondrial genes were removed from the data. This resulted in a total of 9180 cells. Two clusters of Gnai2 were merged as they had very similar expression profile. Solitary chemosensory cells and endothelial cells were manually assigned a cluster identity by selecting the cells using the CellSelector function of Seurat.

To investigate sexually dimorphic gene expression at cluster level, male and female cells of each cluster were compared individually using FindMarkers function of Seurat. For neuron specific analysis, neuronal cell types representing Gnao1 neurons, Gnai2 neurons, immature neurons, and progenitor cells were separated from main seurat object using subset function to create a new ‘neurons’ object. Clustering was performed again using FindClusters function with resolution=0.4 to define neuronal subtypes. This ‘neurons’ object was used for pseudotime analysis with the SlingShot tool. To identify developmental trajectory using SlingShot, progenitor cell cluster(n5) was used as the starting cluster as it expressed Neurod1, Neurog1.

To compare Gnao1 and Gnai2 neurons (mature and immature), differential expression analysis was done using FindMarkers function of Seurat with default parameters. For gene ontology enrichment analysis, differentially expressed genes with log2 fold change greater than 1 were subjected to enrichment analysis using ClusterProfiler’s enricher function. The graphical user interface for scvnoexplorer.com was made using ShinyCell tool (Ouyang et al., 2021).

Co-expression analysis

Non-zero expression level for a particular gene may not indicate that the VR is expressed. Therefore, we identified the cut-off for V1R, V2R and H2-Mvs based on the distribution of expression level across all cells. Normalized gene expression values of V1R, V2R and H2-Mv genes were extracted from each cell of Gnai2 and Gnao1 clusters and the distribution is plotted across all cells. This resulted in a Bi-modal distribution and starting of the second peak was considered as cut-off for V2R/V1R/H2-Mv genes (Figure 7-figure supplement-1D-1F). The genes were considered co-expressed in a single cell if their normalized expression value is greater than the identified threshold value of 1.25 for V2R, H2-Mv gene and 2.5 for V1R genes. Further, the number of cells crossing the threshold were calculated for each combination.

RNA in situ hybridization

Probe design and synthesis

Primers (Supplementary table 9) targeting unique regions of each gene were designed by adding T7/SP6 promoter sequence. PCR (GoTaq DNA polymerase) was performed using gene specific primers with VNO cDNA as template and the product was run on agarose gel to confirm specific amplification of each gene. When multiple bands were seen, the band with desired molecular weight was cut from the gel a. PCR product clean-up or gel purification was performed using NucleoSpin Gel and PCR Clean-up kit. Purified PCR product was verified by sanger sequencing and was used as template for in vitro transcription to obtain digoxigenin (Dig) and fluorescein (Flu) labeled riboprobes using following reaction composition: 1ug Template DNA, 1x Dig/Flu RNA labeling mix, 1U/uL RNAsin plus, 5mM DTT, 1x Transcription buffer in nuclease-free water. The reaction was purified using Qiagen/MN RNA cleanup kit and stored at -80℃ by adding formamide to 50% of volume.

H2-Mv genes and Gnai2 were cloned to a plasmid vector with T7 or SP6 promoters. Gnao1 plasmid was purchased from Invitrogen (6309166). The plasmid was linearized and in vitro transcription was performed as described earlier.

Chromogenic ISH

Fresh VNO was embedded in OCT and 14µm thick sections were collected on a cryostat. Tissue sections were fixed with 4% PFA in Diethyl Pyrocarbonate (DEPC) treated phosphate buffer saline (PBS) and acetylation was performed using Propionic Anhydride, Triethanolamine, NaCl. Permeabilization was done using 0.1M HCl and sections were pre-hybridized using hybridization buffer (50% Formamide, 5X SSC, 5X Denhardt’s solution, 0.1mg/mL Salmon sperm DNA, 0.25mg/mL yeast t-RNA) for 2 hours. Probes were diluted (1:100) in the hybridization buffer and added to the slides. Hybridization was performed for 12-16 hours at 67°C. To remove the unbound probe, sections were washed thrice with 0.2X SSC for 30 minutes each. Slides were equilibrated in a buffer consisting of 0.1M Tris-HCl. pH 7.5, 150mM NaCl for 5 minutes and blocked with 10% FBS in the same buffer. Hybridized Digoxigenin (DIG) or Fluorescein (Flu) containing RNA probes were detected by alkaline phosphatase conjugated Anti-DIG or Anti-Flu Fab₂ fragments (1:7500 dilution). Unbound antibody was washed, and development of alkaline phosphatase was done using BCIP/NBT substrate diluted as per manufacturer’s protocol in 0.1M Tris-HCl pH-9, 0.1M NaCl, 50mM MgCl₂. The signal was developed for 12 - 72 hours based on the intensity. After development, slides were washed with PBS and mounted using 10% glycerol in 0.1M Tris-HCl (pH 7.5).

Two-color ISH

As described above, the same protocol for chromogenic ISH was followed until 0.2X SSC washes, after which the hybridized probes were detected using peroxidase conjugated anti-FLU or anti-DIG Fab(2) fragments and the Tyramide signal amplification (TSA) system from Akoya Biosciences. Slides were incubated in 3% H₂O₂ in PBS for 1 hour and washed thrice for 10 minutes each and blocked with 0.5% blocking buffer for 30 min. Dig and Flu were sequentially developed using TSA-FITC or TSA-Cy3.

Immunohistochemistry

VNOs were dissected from animals of age 8-12 weeks, fixed with 4% Paraformaldehyde (PFA) in PBS and cryopreserved with 30% sucrose. Tissue was embedded in OCT-freezing medium and cryosections of 14µm thickness were collected on glass slides. Sections were post fixed again with 4% PFA in PBS, blocked and permeabilized by incubating for 2 hours with a blocking buffer (3% Bovine Serum Albumin, 0.1% TritonX-100 in PBS with 0.02% Sodium Azide). After blocking, the sections were incubated with primary antibodies diluted in the blocking buffer for 2 hours followed by washing off excess antibody and secondary antibody incubation for 2 more hours. All washes were done using 0.1% TritonX-100 in PBS.

Light microscopy image acquisition

Chromogenic images were acquired on Olympus BX43 upright microscope using bright field Koehler illumination and 10X objective (PlaN, 0.24 NA) equipped with Olympus DP25 color CCD camera. Two color fluorescent RNA-ISH images and fluorescent Immunohistochemistry images were acquired using Leica Stellaris upright confocal microscope using 10X (HC PL APO 0.40 NA) or 63X (HC PL APO 1.40-0.60 NA Oil) objective keeping 1 Airy unit as pinhole size.

Electron Microscopy

VNOs were dissected from animals that were trans-cardially perfused with PBS followed by 4% PFA in PBS buffer. Dissected VNO was embedded in 2% agarose and coronal sections of 400um were obtained on a vibratome in 0.1M sodium cacodylate buffer (pH 7.4). Sections were fixed with Karnovsky’s fixative (3% PFA + 2% glutaraldehyde in 0.1 M sodium cacodylate) for up to 1 week. Subsequent processing steps were similar to that described before (Terasaki, Brunson and Sardi, 2020). Briefly, vibratome sections were rinsed with 0.1 M sodium cacodylate then immersed in 1% osmium, 0.8 % potassium ferricyanide in the cacodylate buffer for 1 hour. This was followed by incubation in 1% aqueous uranyl acetate 1 hr, then 30 min in lead aspartate (Walton, 1979), dehydration in graded ethanol, infiltration with epon resin then embedding and polymerization at 60 °C for 2 days. Serial 70 nm thick sections were collected using a Powertome automated tape collector (RMC Boeckeler, Tucson, AZ). The tape was mounted on a silicon wafer, carbon coated and imaged with a field emission electron microscope (Thermo Fisher / FEI Verios, Hillsboro, OR). Backscatter electrons were imaged from a 5 keV / 0.8 nanoAmp beam of electrons.

Materials

Software

Data availability

All raw data related to single cell RNA sequencing was deposited to GEO and can be publicly accessed using accession ID GSE253252.

Acknowledgements

We acknowledge Jyoti Rohilla, Nandana Nanda for help in preparation of DNA templates for riboprobe generation, Tulasi Nagabandi for helping with scRNA library preparation, Tamal Das for sharing Atlastin1, Sec61β antibodies and Erik Snapp for insightful discussions. We acknowledge support from the Department of Atomic Energy, Government of India, under Project No. RTI 4007.