Hierarchical deconstruction of mouse olfactory sensory neurons: from whole mucosa to single-cell RNA-seq

Abstract

The mouse olfactory mucosa is a complex chemosensory tissue composed of multiple cell types, neuronal and non-neuronal. We have here applied RNA-seq hierarchically, in three steps of decreasing cellular heterogeneity: starting with crude tissue samples dissected from the nose, proceeding to flow-cytometrically sorted pools of mature olfactory sensory neurons (OSNs), and finally arriving at single mature OSNs. We show that 98.9% of intact olfactory receptor (OR) genes are expressed in mature OSNs. We uncover a hitherto unknown bipartition among mature OSNs. We find that 19 of 21 single mature OSNs each express a single intact OR gene abundantly, consistent with the one neuron-one receptor rule. For the 9 single OSNs where the two alleles of the abundantly expressed OR gene exhibit single-nucleotide polymorphisms, we demonstrate that monoallelic expression of the abundantly expressed OR gene is extremely tight. The remaining two single mature OSNs lack OR gene expression but express Trpc2 and Gucy1b2. We establish these two cells as a neuronal cell type that is fundamentally distinct from canonical, OR-expressing OSNs and that is defined by the differential, higher expression of 55 genes. We propose this tiered experimental approach as a paradigm to unravel gene expression in other cellularly heterogeneous systems.

Figure 1. Differential expression analysis of mouse olfactory sensory neurons (OSNs) and whole olfactory mucosa (WOM).

(A) Schematic of the RNA-seq experimental strategy. After dissection of the WOM of OMP-GFP (+/−) male and female mice, pools of ~10 million OSNs were collected by FACS. RNA was extracted from these, along with WOM samples, cDNA generated, and libraries were amplified for deep sequencing. (B) Differential gene expression analysis between the transcriptomes of OSNs and WOM. Statistically significant differentially expressed genes (fold-change >3; FDR < 5%) are highlighted in pink. (C) Comparison of OR and TAAR gene expression levels. A scatter plot of OR and TAAR gene expression levels (black) in the WOM versus the sorted OSNs reveals a strong correlation. The red line represents the 1:1 diagonal. Classical marker genes for mature OSNs (yellow) are similarly enriched. (D) Distribution of OR and TAAR gene expression (normalized counts) as a barplot, for both OSNs and WOM. Genes are displayed in ascending order of their expression values in OSNs (dark green). Matching values are plotted for WOM (bright green), demonstrating that the dynamic distribution of OR gene expression is largely conserved, but differs in OR genes expressed at low levels in WOM (arrow).

Figure 2. Differential expression analysis of mature OSN subpopulations.

(A) Schematic of the RNA-seq experimental strategy. The WOM of OMP-GFP (+/−) male mice was dissected and pools of ~10,000 OSNs from two non-overlapping subpopulations, OMP-GFP^low (light green) and OMP-GFP^high (dark green), were collected by FACS gating based on fluorescence intensity. Subsequently the RNA was extracted, cDNA generated, and libraries were amplified for deep-sequencing. (B) Representative coverage plot showing the distribution of RNA-seq reads across the Omp locus (which is located within an intron of Capn5) is similar in both subpopulations. (C) The number of Omp normalized counts is 1.5 times higher in the OMP-GFP^high population (unpaired t-test, 2 tails, *P = 0.025). Analysis of the same samples by TaqMan qRT-PCR verified this increase. (D) Scatter plot showing the expression levels of previously reported immature (brown) and mature (green) OSN-specific markers in the OMP-GFP^low and OMP-GFP^high OSNs. Both sets of genes are expressed at equivalent levels in both populations (Wilcoxon rank sum test, P = 0.45 and P = 0.22 for mature and immature gene sets), but mature specific genes are expressed considerably higher (Wilcoxon rank sum test, P < 2.2–16 and P < 2.2–16 for OMP-GFP^low and OMP-GFP^high). (E) Differential gene expression analysis between the transcriptomes of OMP-GFP^low and OMP-GFP^high OSNs. Statistically significant differentially expressed (FDR < 5%) genes are highlighted in pink. (F) Functional terms enrichment analysis between the OMP-GFP^low and OMP-GFP^high differentially expressed genes identified a total of 88 GO terms significantly overrepresented (FDR < 5%). The bar graph indicates the expected and observed number of DE genes for six of these GO terms (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001). The full list and GO category names is in Supplementary Data S2. (G) Projection of the OSNs, OMP-GFP^low and OMP-GFP^high samples onto the cell cycle Principal Component Analysis (PCA) that differentiates between cell cycle stages (see Supplementary Fig. S2B and Supplementary Methods for details). Principal Component 1 (PC1) segregates the samples based on their cell cycle stage, with samples in G1 having negative values and samples in S/G2-M positive values. All tested samples are in the G1 stage.

Figure 3. RNA-seq in single mature OSNs.

(A) Schematic of the single-cell RNA-seq experimental strategy. After dissection of WOM of a heterozygous OMP-GFP mouse, 5,000 mature OMP-GFP^high OSNs were collected by FACS into a C1 IFC (10–17 μm, Fluidigm). Using the Fluidigm C1 platform, single cells were sorted, imaged, RNA extracted, cDNA generated, and libraries amplified for deep-sequencing. A representative bright field and fluorescence image for a captured single OSN is shown in the middle panel. (B) Capture rate for single cells: 30 wells captured more than one cell and/or debris, 58 captured single cells, and 8 were empty. Representative bright field images are shown for each of these categories. (C) The 58 captured single OSNs were subjected to stringent, hierarchical quality control criteria, which defined low-quality (LQ, black dots) and high-quality cells (HQ, green dots, see also Supplementary Fig. S3; and Supplementary Methods). (D) Spearman correlation coefficients of the 21 single cells used in our downstream analysis with the WOM and OSN populations (10 million, OMP-GFP^low and OMP-GFP^high OSNs). As expected, the single cells correlate best with the OMP-GFP^high OSNs, followed by the OMP-GFP^low OSNs, the 10 million OSNs and finally WOM. To the right of the graph, boxplots show the cumulative data for each comparison. The thick black horizontal bar corresponds to the median and data points that are outside 1.5 times the interquartile range (the box) are indicated as outliers (black dots).

Figure 4. Captured single cells express canonical marker genes of mature OSNs.

Heatmap of the expression of canonical olfactory cell-specific markers across the WOM and all OSN samples analyzed. While the WOM and OSN population samples express markers for all cell types, the single cells specifically express high levels of mature OSN markers. mOSNs: mature OSNs, iOSNs: immature OSNs, GBCs: globose basal cells, HBCs: horizontal basal cells, SUSs: sustentacular cells, VSNs: vomeronasal sensory neurons.

Figure 5. OR gene expression patterns single OSNs.

(A) Heatmap of the expression levels of the most variable genes in single OSNs. All 598 highly expressed genes (≥1000 normalized counts in at least 1 sample) with a coefficient of variation ≥4 are represented. The identity of these genes are presented in Supplementary Data S3. (B) Normal-like distributions were fitted for all lowly expressed (blue line) and highly expressed (red line) OR genes. Only genes at or above the intersection of both curves (855 normalized counts) were considered as abundantly expressed. (C) Expression level of the most abundant OR in each of the 21 single OSNs. The most abundant ORs in OSN 259 and OSN 261 are expressed at extremely low levels. (D) Phylogenetic tree of all mouse OR genes (Class I and II in black and grey, respectively). The most abundant ORs in all the single OSNs analyzed are depicted as red circles. (E) Expression levels of all ORs in OSN 188. The two most abundant ORs in that cell are indicated. (F,G) Coverage plots of the RNA-seq data for OSN 188 mapping to the first (Olfr1348, (F)) and second (Olfr1372-ps1, (G)) most abundant ORs. Boxes correspond to exons and arrowheads indicate the strand of the gene. The existing Ensembl annotations are shown in the top box. Optimized OR gene models used in our analysis are shown in the middle box. The sequencing data mapping to the OR gene models is below. (H) Expression levels of all ORs in OSN 236. The two most abundant ORs in that cell are indicated. (I,J) Coverage plots of the RNA-seq data for OSN 236 mapping to the first (Olfr556, (I)) and second (Olfr1509 (J)) most abundant ORs. Boxes correspond to exons and arrowheads indicate the strand of the gene. The existing Ensembl annotations are shown in the top box. Optimized OR gene models used in our analysis are shown in the middle box. The sequencing data mapping to the OR gene models is below.

Figure 6. Monoallelic expression of OR genes in single OSNs.

(A,B) Coverage plots of the sequencing data mapping to the same OR (Olfr55) in two cells, OSN 222 (A) and OSN 263 (B). Boxes correspond to exons and arrowheads indicate the strand of the gene. The existing Ensembl annotations are shown in the top box. Optimized OR gene models used in our analysis are shown in the middle box. The sequencing data mapping to the OR gene models is below. Reported single nucleotide polymorphisms (SNPs) used to map the C57BL/6 and 129P2 alleles are shown in black and golden vertical lines, respectively. (C,D) The percentage (%) of sequencing data showing the C57BL/6 (black) or the 129P2 (gold) allele at each SNP on Olfr55 in OSN 222 (C) and OSN 263 (D). Genomic coordinates on chromosome 17 for each SNP are indicated in the x-axis. (E) The percentage (%) of sequencing data showing the C57BL/6 (black) or the 129P2 (gold) allele for all ORs with at least one reported informative SNP. In all cases, only one allele is expressed per cell. (F) Example of differential splicing of an OR gene. The sequencing reads in the bottom panel are represented by grey boxes; blue lines join portions of reads that map across exon junctions. Two isoforms are present, with the differential inclusion of the second exon. Thus individual OSNs express multiple isoforms of the chosen OR gene.

Figure 7. Gene expression profile of another neuronal type in the MOE.

(A) Scatter plot of the mean expression values for all genes in the two cells lacking OR gene expression versus the remaining 19 OSNs. Statistically significant differentially expressed genes (FDR < 5%) are highlighted in pink; those that show consistent high expression in both cells lacking OR gene expression are in purple. The three most abundant DE genes are indicated (Gucy1b2, Sln and Emx1). (B) Heatmap of the 55 DE genes that constitute the gene expression signature of the cells lacking OR gene expression. (C,D) Cryosections of adult mouse MOE hybridized with cRNA probes for the two most highly expressed DE genes – Gucy1b2 (C) and Sln (D). The hybridization signals are sparsely distributed within the MOE. (E) Two-color in situ hybridization of the top ranked marker (Gucy1b2) with Trpc2. As previously shown, some, but not all of the Trpc2 cells in the MOE also coexpress Gucy1b2. (F) Two-color in situ hybridization of the top ranked marker (Gucy1b2) with other differentially expressed genes, Sln, Emx1 and Sncg. Arrowheads point to labeled cells. Scale bars, 50 μm.