Single cell transcriptome analysis of mouse carotid body glomus cells

Abstract

Key points: Carotid body (CB) glomus cells mediate acute oxygen sensing and the initiation of the hypoxic ventilatory response, yet the gene expression profile of these cells is not available. We demonstrate that the single cell RNA-Seq method is a powerful tool for identifying highly expressed genes in CB glomus cells. Our single cell RNA-Seq results characterized novel CB glomus cell genes, including members of the G protein-coupled receptor signalling pathway, ion channels and atypical mitochondrial electron transport chain subunits. A heterologous cell-based screening identified acetate (which is known to affect CB glomus cell activity) as an agonist for the most highly abundant G protein-coupled receptor (Olfr78) in CB glomus cells. These data established the first transcriptome profile of CB glomus cells, highlighting genes with potential implications in CB chemosensory function.

Abstract: The carotid body (CB) is a major arterial chemoreceptor containing glomus cells whose activities are regulated by changes in arterial blood content, including oxygen. Despite significant advancements in the characterization of their physiological properties, our understanding of the underlying molecular machinery and signalling pathway in CB glomus cells is still limited. To overcome this, we employed the single cell RNA-Seq method by performing next-generation sequencing on single glomus cell-derived cDNAs to eliminate contamination of genes derived from other cell types present in the CB. Using this method, we identified a set of genes abundantly expressed in glomus cells, which contained novel glomus cell-specific genes. Transcriptome and subsequent in situ hybridization and immunohistochemistry analyses identified abundant G protein-coupled receptor signalling pathway components and various types of ion channels, as well as members of the hypoxia-inducible factors pathway. A short-chain fatty acid olfactory receptor Olfr78, recently implicated in CB function, was the most abundant G protein-coupled receptor. Two atypical mitochondrial electron transport chain subunits (Ndufa4l2 and Cox4i2) were among the most specifically expressed genes in CB glomus cells, highlighting their potential roles in mitochondria-mediated oxygen sensing. The wealth of information provided by the present study offers a valuable foundation for identifying molecules functioning in the CB.

Figure 1. Workflow of preparing single cell RNA‐Seq

The tissue of interest was dissected out from mice for further enzymatic and mechanical treatment to obtain single cells in suspension. Cells are coloured differently to represent the heterogeneity of cell types present at this step. Each individual single cell was carefully dispensed into a PCR tube containing lysis buffer and heated for RNA denaturation. mRNAs were reverse transcribed using anchor T primers and further processed to generate poly(A) tailed single‐stranded cDNAs, which were subsequently amplified using anchor T primers. To identify the cell type of interest, each single cell‐derived cDNA was used as a template for the marker genes analysis. cDNAs capable of amplifying specific cell type markers were selected and used for Illumina library preparation and sequencing. Sequences from these candidate cells were then aligned to the current mouse genome database for further analysis.

Figure 2. Transcriptome comparison of single OSNs

A, scatter plot of gene expression profiles (RPM values) from individually picked and processed single OSNs. The Pearson correlation coefficient between the two samples is indicated. Coloured dots represent the predominant olfactory receptor expressed in each OSN. B, scatter plot comparing the average expression profile of our single OSNs with the transcriptome of fluorescence‐activated cell sorted mature OSNs (Omp positive cells). The Pearson correlation coefficient is indicated. Coloured dots label genes crucial to OSN function. Omp, olfactory marker protein; Gnal, guanine nucleotide binding protein, α stimulating, olfactory type; Cnga2, cyclic nucleotide gated channel α2; Rtp1, receptor transport protein 1; Adcy3, adenylate cyclase type 3.

Figure 3. Transcriptome comparison of single CB glomus cells

A, representative marker genes analysis for four single cells isolated from CB dissociations. Cells 1 and 2 are considered candidate CB glomus cells as a result of their expression of glomus cell markers Th (tyrosine hydroxylase), Uchl1 (ubiquitin carboxy‐terminal hydrolase L1) and Kcnk3 (potassium channel subfamily K member 3). B, histogram of the RPM values of all detected genes in eight individually picked and processed candidate CB glomus cells (G1 through G8). C, scatter plot of gene expression profiles (RPM values) from two single CB glomus cells. The Pearson correlation coefficient between the two samples is indicated. D, scatter plot of mean expression profiles from single CB glomus cell RNA‐Seq and C57Bl6/J CB RNA‐Seq. The Pearson correlation coefficient is indicated. E, scatter plot of mean expression profiles from single CB glomus cell RNA‐Seq and DBA/2J CB microarray. Mean expression levels were calculated and their relative rankings were used for comparison. The Pearson correlation coefficient is indicated.

Figure 4. Validation of single CB glomus cell RNA‐Seq data using in situ hybridization and immunohistochemistry

A–E, genes identified through single CB glomus cell RNA‐Seq were confirmed by in situ hybridization. CB sections from mice 3 weeks or older were hybridized with DIG‐labelled antisense RNA probes of GPCR signalling components (A), ion channels and associated proteins (B), HIF pathway components (C), neuronal markers (D) and others (E). Olfr588, not detected in CB glomus cell RNA‐Seq, served as a negative control (A). Scale bar = 100 μm. F, representative images of adult mouse CB sections stained with primary antibodies specific for the proteins of several highly abundant genes identified in the single CB glomus cell RNA‐Seq. Scale bar = 100 μm. Gnas, (guanine nucleotide binding protein, α stimulating) complex locus; Rgs5, regulator of G‐protein signalling 5; Dgkh, diacylglycerol kinase eta; Adora2a, adenosine receptor A2a; Rgs4, regulator of G‐protein signalling 4; Olfr78, olfactory receptor 78; P2ry12, P2Y purinoceptor 12; Pkib, cAMP‐dependent protein kinase inhibitor β; Dgkk, diacylglycerol kinase kappa; Lphn1, latrophilin‐1 precursor; Adcy1, adenylate cyclase type 1; Adcyap1r1, pituitary adenylate cyclase‐activating polypeptide type I receptor; Ednra, endothelin‐1 receptor precursor; Hcrtr1, orexin receptor type 1; Olfr558, olfactory receptor 558; Grina, glutamate receptor, ionotropic, N‐methyl‐d‐aspartate‐associated protein 1; Chrna3, neuronal ACh receptor subunit α‐3 precursor; Gria3, glutamate receptor 3; Trpm7, transient receptor potential cation channel subfamily M member 7; Cacnb3, voltage‐gated L‐type calcium channel subunit β‐3; Gria2, glutamate receptor 2; Ndufa4l2, NADH dehydrogenase [ubiquinone] 1α subcomplex subunit 4‐like 2; Epas1, endothelial PAS domain‐containing protein 1; Cox4i2, cytochrome c oxidase subunit 4 isoform 2; Egln1, egl‐9 family hypoxia‐inducible factor 1; Arnt2, aryl hydrocarbon receptor nuclear translocator 2; Hif1α, hypoxia‐inducible factor 1‐α; Nnat, neuronatin; Chga, chromogranin a; Chgb, chromogranin b; Th, tyrosine hydroxylase; Uchl1, ubiquitin carboxy‐terminal hydrolase L1; Cyb561, cytochrome b561; Maged1, melanoma‐associated antigen D1; Methig1, methyltransferase hypoxia inducible domain containing 1; Ly6h, lymphocyte antigen 6 complex, locus H; Adipor1, adiponectin receptor 1; Slc25a4, solute carrier family 25 (mitochondrial carrier, adenine nucleotide translocator), member 4; Car2, carbonic anhydrase 2; Cox4i1, cytochrome c oxidase subunit 4 isoform 1; Npr2, atrial natriuretic peptide receptor 2 precursor; Syp, synaptophysin.

Figure 5. Transcriptome comparison between single CB glomus cells and other cell types

For each cell type, the average RPM values are the mean expression levels from separately picked and processed single cells. A, scatter plot showing a weak correlation between CB glomus cells and OSNs. The Pearson correlation coefficient is indicated. B, scatter plot showing a weak correlation between CB glomus cells and VSNs. The Pearson correlation coefficient is indicated. C, principal component analysis of gene expression patterns for eight single CB glomus cells, two single OSNs, and two single VSNs. The amount of variance explained by principal components 1 and 2 is shown on the x‐ and y‐axes, respectively.

Figure 6. Differential gene expression analysis between single CB glomus cells and non‐CB tissues

Differential gene expression analysis was performed using single cell RNA‐Seq data from CB glomus cells, OSNs, VSNs as well as public available RNA‐Seq data from 15 other mouse tissues. A, volcano plot displays −log₁₀ (P value) as a function of log₂ (fold change) between CB glomus cells and non‐CB tissues, with coloured dots indicating significantly differentially expressed genes (FDR < 0.01). The significantly over‐represented genes in CB glomus cells are indicated in red, whereas the significantly under‐represented genes in CB glomus cells are indicated in blue. B, heat maps showing the expression patterns of the significantly differentially expressed CB glomus cell genes, with the top 50 over‐represented genes shown on the left and top 50 under‐represented genes shown on the right. Mean expression levels were calculated and their relative rankings in each tissue were shown, with the least abundant gene and other non‐detected genes ranked as 100%.

Figure 7. Olfr78 is a conserved olfactory receptor activated by acetate

A, unrooted phylogenetic tree of Olfr78 family orthologs based on similarity of amino acid properties. B, lacZ staining of a carotid artery bifurcation from a heterozygous Olfr78^tm1Mom/MomJ mouse. Blue precipitation marks the lacZ positive cells. SCG, superior cervical ganglion. C, dose–response curves of Olfr78 against SFCA. The x‐axis represents molar concentration in log scale. D, dose–response curves of Olfr78 against acetic acid and sodium acetate (pH 7.4). The x‐axis represents molar concentration in log scale.