Single cell RNA sequencing of the mouse vomeronasal neuroepithelium

A) Uniform Manifold Approximation and Projection (UMAP) of 9185 cells from vomeronasal neuroepithelium with 18 identified clusters. Each point represents a cell that is color coded according to the cell type. Clusters were assigned a cell type identity based on previously known markers. Percentage of total cells corresponding to each cluster is in Supplementary table 3. B) Dot plot showing the expression of top 10 gene markers for each cluster identified by differential expression analysis. Gene expresion markers shared across multiple clusters are listed only once. The dot radius and dot color indicate the percentage of cells expressing the gene and scaled expression value in that cluster. respectively. Top 20 gene markers for each cluster with log2 fold change are listed in Supplementary table 2. C) Feature plot showing the expression level of major known gene markers of neuronal and sustentacular cell types.

Gene markers of non-sensory cells in the vomeronasal organ identified from single cell RNA sequencing.

A) UMAP feature plot (top) of a single representative gene marker for non-sensory cell clusters and their spatial gene expression pattern by RNA in situ hybridization on VNO coronal sections (bottom, scale bar: 100µm). B, C) Three custers of macrophage-like cell types represented by UMAP projection (B) and a dot plot (C) showing differ­ entially expressed genes between these clusters. The dot radius indicates percentage of cells expressing a gene, whose scaled expression level is represented by intensi­ ty of color as per the scale. Additional markers of each non-sensory cell type are in Supplementary table 1 and differential expression between macrophage clusters is in Supplementary table 4.

Gene expression dynamics during development of sensory neurons

A) UMAP of neuronal cells with clusters (n1-n13 from Figure 3, Figure 3-figure supplement-1) represented in different colors. Sub-clusters within mature Gnai2, Gnao1 neurons are merged_ The line on the UMAP plot shows pseudotime developmental trajectory of neurons. BJ Volcano pot of differential gene expression between immature Gnai2 (cluster n6; Gap43+,Gnai2+) vs immature Gnao1 (cluster n7; Gap43+,GnaO1+) neurons_ Genes that satisfy 11092 fold change normalized expression! > 1 (green) and -10910 p value > 6 (blue) are considered significant and coloured in red_ Non-significant (NS) genes are colored in grey_ Arrows point to transcription regulators enriched in Gnao1+ immature neurons. Complete list of differentially expression genes is in Supplementary table 6. CJ Feature plots showing the normalized expression levels of previously known markers for immature neurons (Gap43), Gnao1 neurons (Tfap2e), Gnai2 neurons (Meis2) and transcription factor or related genes: Creb5, Prrxl1, Shisa8, Lmo4, Foxo1_Arrows highlight the limited expression of Creb5 and Prrxl1 to immature neurons, but absent from mature indicating transient expression during development of Gnao1 neurons. Subclusters within Gnao1, Gnai2 neurons and effect of VRs on subclustering are in Figure 3-figure supplement 1, 2 respectively.

Subpopulation of Gnao1 neurons defined by V2R and H2-Mv expression

A) UMAP representation of Gnao1 neurons from Figure 3. Each dot represents a cell and four Gnao1 neuron clusters (n1-n4) are marked in different colours. B-D) Feature plot showing exclusive expression of Vmn2r1 (B), Vmn2r2 (CJ and the most abundant H2-Mv gene, H2-M10.3 (DJ in Gnao1 neurons. E) Heat map showing the expression of V2R and H2-Mv genes in Gnao1 neurons. Cluster numbers are marked on the top with color coding as in (A) and gene families are labelled on the left. Each column represents a cell and the scaled gene expression in each row is colour coded as per the scale with red and blue indicating high and low expression respectively. Vmn2r1 is expressed in almost all cells of cluster-1 and cluster-4; Vmn2r2 is expressed in all cells of cluster 3; Cluster2 has mutually exclusive expression of Vmn2r1 and Vmn2r2. F-H) Bar plot showing number of cells expressing: 0-7 family-C V2Rs per cell (F), 0-5 family-ABD V2Rs per cell (G), 0-6 H2-Mv genes per cell (H) with composition of cells associated with family C1 (orange) or C2 (blue) V2R color coded on the bar. I) Box plots comparing the sum of normalised expression levels of family-C V2Rs and Gnao1 (Green) in a cell that expresses 1 to 5 fami­ ly-C V2Rs. J) Box plot comparing the level of total ABD-V2R expression from cells with 1 and 2 ABD-VRs along with Gnao1 expression level (green). K) Box plot comparing the level of total H2-Mv expression in cell that express 1-5 H2-Mv genes along with Gnao1 expression level (Green). Multiple combinations of family-C, family ABD V2Rs and H2-Mvs identified to be co-expressed in a single cell and their cell frequency are listed in Supplementary table 7.

Divergent gene expression: Gnao1 vs Gnai2 neurons.

A) Volcano plot showing differentially expressed genes between mature Gnao1 and Gnai2 neurons. Genes that satisfy |log2 fold change normalized expression| > 2 (green) and -log10 p value > 6 (blue) are considered significant and coloured in red. Non-signifi- cant (NS) genes are colored in grey. Complete list of differentially expressed genes is in Supplementary table 6. C-D)Two color RNA ISH of selected enriched genes marked in bold on the volcano plot. Gnao1, Gnai2, respective markers of basal and apical neurons are shown in (B). Genes enriched in Gnai2 neurons (C) or Gnao1 neu- rons (D) are co-localized with the respective markers. RNA-ISH for additional enriched genes are shown in Figure 5-figure supplement-1. Scale bar:100 μm.

A) ER gene expression in Gnao1 neurons.

Annotation of gene ontology (GO) biological processes of genes that are significantly enriched in Gnao1 neurons from Figure 5A. GO terms related to ER processes are marked in red. p-value < 0.05 was used as cut-off. B) VNO coronal section two color fluorescent RNA-ISH of selected Gnao1 enriched ER chaperone genes (Creld2, Pdia6, Dnajc3, Sdf2I1: green), co-labelled with Gnao1 probe (red) shows Gnao1 zone restricted expression of these genes. Scale bar 100µm.

Differential localization of ER proteins in VNO neurons.

Pseudocolored immunofluorescence images of VNO coronal sections labelled with antibodies against KDEL (anti-SEKDEL), a common ER retention signal (A). The section is co-labelled with anti-Gnao1 to mark basal zone neurons (B) and anti-OMP to mark all mature neurons (C). Signal intensity of KDEL, Gnao1, OMP channels is quantified from ROIs along the apical-basal axis as shown in example (D). Signal intensity measured from multiple sections (n>20 for each anitbody) from 3 animal replicates was normalized and trendline was fitted to show the Gnao1 neuron enriched localization of anti-KDEL (E). Points on the plot show normalized intensity values color coded for each antibody on which the trendline was fitted. Similar imuno-labelling and quantification of ER chaperone Hspa5 (BiP) (F), ER membrane WUDQVRORFRQ VXEXQLW 6HF61ȕ (G) and ER membrane protein Atlastin1 (H) indicate their enrichement in Gnao1 neurons compared to Gnai2 neurons. Distribution of additional ER chaperone and membrane proteins is shown in Figure 7-figure supplement-1. Scale bar: 50 μm.

Basal/Gnao1 neurons are densely packed with cubic membrane ER.

A) Scanning electron micrographs of vomeronasal sensory epithelium at low magnification showing cell bodies of VSNs, sustentacular cells and basal lamina. Boxed regions in red or green mark cell bodies of basal/Gnao1 or apical/Gnai2 neurons respectively, that are displayed at higher magnifica- tion below. Nucleus (N) appears light and ER is dark. Cell bodies of basal/Gnao1 neurons are larger and occupied by substantial amount of ER, incomparison to apical/Gnai2 neurons. B, B’) Magnified micrographs show the cell body of a basal/Gnao1 neuron, packed with cubic ER membranes. White arrows point to dense membrane stacks that are better resolved in Figure 8-figure supplement 1. Red arrow points to lamellar ER membrane that is contiguous with the cubic membrane. C, C’) Cell bodies of apical/Gnai2 neurons also seem to show dense cubic membrane ER, that is smaller in comparison to basal/Gnao1 neurons.

Onset of V2R expression coincides with the expression of ER chaperone genes

A) Feature Plot showing the expression of Gnao1, Vmn2r1, Sdf1I1 and Mani. Sdf2I1 and Mani are known ER chaperones and their upregulation in Gnao1 neurons coincides with Vmn2r1 expression, which is preceded by Gnao1 expression. B) Heatmap showing the expression of Gnao1, Vmn2r1, Vmn2r2 and sev­ eral ER chaperone genes (red) in the clusters arranged as per their developmental trajectory. C) Cartoon summarizing major transcription factor expression during development leading to Gnao1 neurons with chaprone rich hypertropic ER compared to Gnai2 neurons.

Comparison of cell type composition and gene expression from male and female vomeronasal neuroepithelium.

A) Uniform Manifold Approximation Projection (UMAP) of cells from male and female vomeronasal neuroepithelium. with the cluster numbers corresponding to Figure 1A, 1B. Solitiary chemosensory neuron (Cluster 18) were seen only in female data. B) Scatter plots comparing average expression of genes accross each cluster from male and female with Pearson correlation co-efficient at the top of the plot. Each point in the plot represents a gene. Known sexually dimorphic genes: Eif2s3y. Ddx3y, Uty. Kdm5d are marked in red. Scatter plot of cluster 17 between male and female is not shown due to low cell count. The results of differential expression between each cluster of male and female are in Supplementary table 1.

A) UMAP projection of neurons with 13 clusters (n1-n13). B) Feature plot showing expression of neuronal markers associated with various stages of differentiation: Globose basal cells (n5; Ascl1), progenitors cells (n5; Neurod1, Neurog1), immature neurons (n6, n7; Gap43+), mature Gnao1 neurons (n1-n4; Gap43-, Omp+, Gnao1+) and mature Gnai2 neurons (n8-n13; Gap43-, Gnai2+, Omp+). Arrows highlight the expression of Gnao1 and Gnai2 in Gap43+ immature neurons (n6, n7). C) Dotplot showing enriched genes in each cluster compared to all other clusters obtained by differential expression analysis. Size of the dot represents the percentage of cells expressing the gene in that cluster and colour indicates the scaled expression value. Top five gene markers based on log2fold change from each cluster were chosen by filtering the genes whose adjusted p-value is less than 0.005, expression in atleast 50% of cells and less than 50% of cells of all other clusters. No markers were found for cluster n9 with this filtering criteria. Complete list of top 30 markers for neuronal clusters is in Supplementary table 5.

Effect of VR genes on neuronal clustering.

Clustering of neurons based on top 2000 variable geneset in the dataset leading to 13 clusters (same as Figure 3-figure supplement 1A). Reclustering performed after excluding genes coding vomeronasal receptors from the variable geneset without changing any other parameters. The bifurcation of Gnao1, Gnai2 neurons at mature and immature stages remains unchanged. The only changes are merger of mature Gnao1 sub-clusters n1/n3 and change in mature Gnai2 n8 composition.

Bar plot showing the number of cells expressing top 20 V1Rs (A), V2Rs (B) and H2-Mvs (C). Normalized gene expression values were extracted for all V1Rs (D) from Gnai2 neurons, V2Rs (E) and H2-Mvs (F) from Gnao1 neurons and density was plotted to identiy the distribution of gene expression. The red dashed line represents a normalized gene expression value that was used as threshold to call the expression of respective genes.

Characteristics of H2-Mv expression

A -B) Feature plot showing the expression H2-M10 family genes (A) and limited expression of phylogentically divergent H2-Mv members H2-M1, H2-M9 and H2-M11 (B) in Gnao1 neurons. C) Feature plot showing the expression of Vmn2r81 and Vmn2r82 in few cells of cluster 2 and 4 of Gnao1 neurons that express H2-M1, M9 and M11. D) Scatter plot showing the normalized expression level per cell of rarely expressed H2-Mv genes (H2- M1, M9 and M11) on x-axis and Vmn2r1 or Vmn2r2 on y-axis indicating that H2-M1, M9 and M11 co-express with Vmn2r1 unlike other H2-Mv genes. E) Two-color RNA in situ hybridization of H2-M1,H2-M9 and H2-M11 with Vmn2r81/82 confirming the co-expression. The ISH probe was common for 9U81 DQG 82; +2-01/09/011 LQGLFDWHV SRROLQJ RI LQGLYLGXDO SUREHV. 6FDOH EDU: 100 ȝP.

Coexpression characteristics of V2R and H2-Mv genes using data from Hills et al., 2024

A-C) Bar plot showing number of cells expressing: 0-7 family-C V2Rs per cell (A), 0-5 family-ABO V2Rs per cell (8), 0-6 H2-Mv genes per cell (C) with composition of cells associated with family C1 (orange) or C2 (blue) V2R color coded on the bar. The trend of co-expression is similar to Figure 7F, 7G and 7H indicating that co-expression counts are similar across datasets.

Co-expression of V1Rs in Gnai2 neurons

A) Heatmap showing the expression of selected V1Rs in Gnai2 neurons. Each column represents a cell and the scaled expression value of each VR is color coded in each row with red and blue indicating high and low expression respectively. A continues red line in two rows of a single column indicates the expression of two receptors in a single cell. B) Bar plot showing number of cells expressing 0-3 V1Rs per cell indicating the V1R co-expression is limited to small subset of cells. C) Box plot comparing total V1R and Gnai2 (green) expression from cells shown in B. Multiple combinations of V1Rs co-expressed in a single cell and their cell frequency are listed in Supplementary table 8.

Expression pattern of enriched genes in Gnai2 or Gnao1 neurons.

Chromogenic RNA-ISH showing expression of Gnai2 enriched genes: Nsg1, Rtp1, Dner, Qpct, Pcdh7, Prph (A) and Gnao1 enriched genes: Apmap, Selenof, Hspa5 (Bip), Itm2b, Agpat5, Sncg (B). Sncg and Prph are expressed in a scattered pattern amongst few neurons. Scale bar: 100 μm

Pseudocolored immunofluorescence images of VNO cornonal sections labeled with antibodies against ER chaperone proteins: PDI / Grp94 / Calnexin (A), ER structural proteins: NogoB / Climp64 / Reep5 (B). Higher intensity in Gnaol neurons compared to Gnai2 neurons is seen from the images and quantified along the apical-basal axis by co-labelling sections with anti-OMP to mark all neurons and anti-Gnao1 (or anti-SEKDEL -see Figure 7A-E epending on antibody species), to mark Gnaol neurons. Fluorescence intensity along the apical-basal axis is quantified from atleast 20 sections from 3 nimal replicates. The normalized fluorescence intensity of ER proteins increases along apical-basal axis, similar to Gnaol or SEKDEL, while the level of MP either remains same or slightly decreases. Scale bar: 50 pm

Comparison of ER gene expression between Gnai2, Gnao1 neurons.

Violin plots showing the gene expression levels in mature Gnao1 or Gnai2 neurons whose protein levels are shown via immunofluorescence in Figure 7 and Figure 7-figure supplement 1. Log, fold change value for Gnao1 vs Gnai2 calculated from pseudobulk differential gene expression analysis is mentioned on top of the line and Bonferroni-adjusted p-value is mentioned below the line (ns - not significant). RNA levels of Hspa5. Calnexin. Hsp90b1. Sec61b. Reep5 are significantly upregulated in Gnao1 neurons while POI and Atlastin1 do not differ significantly between Gnao1 and Gnai2 neurons.

Scanning electron micrographs showing cell bodies of basal/Gnao1 neurons with sinusoidal ER membranes interspersed with stacked membranes (arrows).