Biology and Bias in Cell Type-Specific RNAseq of Nucleus Accumbens Medium Spiny Neurons

Abstract

Subcellular RNAseq promises to dissect transcriptional dynamics but is not well characterized. Furthermore, FACS may introduce bias but has not been benchmarked genome-wide. Finally, D1 and D2 dopamine receptor-expressing medium spiny neurons (MSNs) of the nucleus accumbens (NAc) are fundamental to neuropsychiatric traits but have only a short list of canonical surface markers. We address these gaps by systematically comparing nuclear-FACS, whole cell-FACS, and RiboTag affinity purification from D1- and D2-MSNs. Using differential expression, variance partitioning, and co-expression, we identify the following trade-offs for each method. RiboTag-seq best distinguishes D1- and D2-MSNs but has the lowest transcriptome coverage. Nuclear-FACS-seq generates the most differentially expressed genes and overlaps significantly with neuropsychiatric genetic risk loci, but un-annotated genes hamper interpretation. Whole cell-FACS is more similar to nuclear-FACS than RiboTag, but captures aspects of both. Using pan-method approaches, we discover that transcriptional regulation is predominant in D1-MSNs, while D2-MSNs tend towards cytosolic regulation. We are also the first to find evidence for moderate sexual dimorphism in these cell types at baseline. As these results are from 49 mice (n_male = 39, n_female = 10), they represent generalizable ground-truths. Together, these results guide RNAseq methods selection, define MSN transcriptomes, highlight neuronal sex differences, and provide a baseline for D1- and D2-MSNs.

Figure 1

Library characterization demonstrates fewer differentially expressed transcripts and a predominance of protein coding genes in RiboTag compared to whole cell and nuclear RNAseq. (a) Method scheme for whole cell FACS, nuclear FACS, and RiboTag affinity purification showing key differences in steps involving sample preparation, cell dissociation, and retrieval of cellular fractions. (b) Density plot for all methods of ln([average FPKM] + 1) of all genes captured shows that RiboTag does not capture a number of genes. (c) Density plots for all methods of ln([average FPKM] + 1) of differentially expressed genes show similar distributions across method. (d) Overlap of D1- and D2-enriched differentially expressed genes across all methods (total D1 overlap = 134; total D2 overlap = 64). (e) Fold change of differentially expressed genes from most D1-enriched in yellow to most D2-enriched in blue, sorted by fold change in whole cell (black = low fold change, grey = not detected in the dataset). (f) Density plots for each method of ln([average FPKM] + 1) of D1 (dark line) and D2 (light line) differentially expressed genes in the respective cell types; medians are indicated with dashed lines. (g) Mean-variance plots comparing ln(variance) to ln([average FPKM] + 1) of differentially expressed genes in pooled D1- and D2-MSNs show only slight differences across methods. (h) Gene biotype distributions for each method’s differentially expressed genes show a decreasing proportion of protein coding genes from the RiboTag to the whole cell to the nuclear dataset.

Figure 2

Method accounts for the most variance and method-variable genes are associated with their expected subcellular compartments. (a) Principal component analysis (PCA) across all samples showing separation by method; PCA within method showing separation by cell type (insets). (b) Percent variance explained by method, cell type and residuals with GO enrichment of method-, cell type-, and residual-variable genes. (c) Gene expression profiles of the top 20 and top 10 most method- and cell type-variable genes, respectively. (d) Subcellular localization of method-variable genes shows enrichment of nuclear-variable genes in the nucleus, and RiboTag-variable genes in the ribosome, as expected.

Figure 3

Method-correlated WGCNA modules provide insight into the biological function of key regulatory and unannotated genes. (a) Dendrogram of hierarchical relationship and cluster-map showing correlation values of module and method eigengenes. (b) Whole cell hub network including top 20 hub genes (labeled) and their first degree edges; enrichment of GWAS genes among whole cell hub network genes (Fisher’s Exact Test). (c) RiboTag hub network including top 20 hub genes (labeled) and their first degree edges; enrichment of GWAS genes among ribosomal hub network genes. (d) Nuclear hub network including top 20 hub genes (labeled) and their first degree edges; enrichment of GWAS genes among nuclear hub network genes.

Figure 4

D1 hub networks from multiple methods demonstrate convergent biological function and transcriptional regulation. (a) WGCNA cluster-map with correlations among module eigengenes and D1 cell type; D1 was most correlated with the cyan module in whole cell (Pearson correlation coefficient = 0.91; p-value = 7.0 × 10⁻⁷), the brown module in nuclear (corr = 0.81; p-value = 2.0 × 10⁻³), and the turquoise module in RiboTag (corr = 0.94; p-value = 8.0 × 10⁻⁶). (b) D1 module hub networks for each method with hub nodes labeled. (c) D1 module GO and KEGG/reactome pathway terms with −log(p-value) of enrichment. (d) Overlap of top 100 whole cell, nuclear and RiboTag TFBSs with number of corresponding transcription factors indicated above the arrow.

Figure 5

D2 hub networks from multiple methods demonstrate convergent biological function and transcriptional regulation. (a) WGCNA cluster-map with correlations among module eigengenes and D2 cell type; D2 was most correlated with the purple module in whole cell (corr = 0.83; p-value = 7.0 × 10⁻⁵), the turquoise module in nuclear (corr = 0.86; p-value = 6.0 × 10⁻⁴), and the blue module in RiboTag (corr = 0.98; p-value = 6.0 × 10⁻⁸). (b) D2 module hub networks for each method with hub nodes labeled. (c) D2 module GO and KEGG/reactome pathway terms with −log(p-value) of enrichment. (d) Overlap of top 100 whole cell, nuclear and RiboTag TFBSs with number of corresponding transcription factors indicated above the arrow.

Figure 6

RiboTag data can discern small-variance biological variables, with WGCNA generating sex-correlated modules. (a) RiboTag WGCNA cluster-map showing correlation values of module eigengenes and sex; male sex was most correlated with the black module (corr = 0.27; p-value = 0.2) and female sex was most correlated with the pink module (corr = 0.55; p-value = 8.0 × 10⁻³). (b) Black module hub network (male) with hub nodes labeled; GO and pathway terms with −log(p-value) of enrichment displayed. (c) Pink module hub network (female) with hub nodes labeled; GO terms and TFBS prediction with −log(p-value) of enrichment displayed. (d) Overlap of D1- and D2-DEGs from male and female datasets; union heatmap of these DEGs showing log₂foldchange (D1 vs. D2) in each dataset (black = low fold change, grey = not detected in the dataset).