TARGET-seq: Linking single-cell transcriptomics of human dopaminergic neurons with their target specificity

Abstract

Dopaminergic (DA) neurons exhibit significant diversity characterized by differences in morphology, anatomical location, axonal projection pattern, and selective vulnerability to disease. More recently, scRNAseq has been used to map DA neuron diversity at the level of gene expression. These studies have revealed a higher than expected molecular diversity in both mouse and human DA neurons. However, whether different molecular expression profiles correlate with specific functions of different DA neurons or with their classical division into mesolimbic (A10) and nigrostriatal (A9) neurons, remains to be determined. To address this, we have developed an approach termed TARGET-seq (Tagging projections by AAV-mediated RetroGrade Enrichment of Transcriptomes) that links the transcriptional profile of the DA neurons with their innervation of specific target structures in the forebrain. Leveraging this technology, we identify molecularly distinct subclusters of human DA neurons with a clear link between transcriptome and axonal target-specificity, offering the possibility to infer neuroanatomical-based classification to molecular identity and target-specific connectivity. We subsequently used this dataset to identify candidate transcription factors along DA developmental trajectories that may control subtype identity, thus providing broad avenues that can be further explored in the design of next-generation A9 and A10 enriched DA-neurons for drug screening or A9 enriched DA cells for clinical stem cell-based therapies.

Keywords: A9 and A10 dopamine neurons; Axonal projections; Parkinson’s disease; human embryonic stem cells; single-cell sequencing.

Fig. 1.

Cell type composition of 12-mo-old fetal and hPSC-DAergic grafts. (A) Amphetamine-induced rotation test indicating functional recovery in transplanted human fetal cell (n = 12) and VM-patterned hPSCs (n = 10) after 6 and 12 mo (one-way ANOVA: ∗∗∗P < 0.0001 for pre-Tx vs. 6 and 12 mo post-Tx). (B) Immunohistochemistry of TH/hNCAM in human fetal cell, and (C) VM-patterned hPSCs intranigral grafts 12 mo posttransplantation. (Scale bar, 200 µM.) (D) Analysis of graft volume 12 mo posttransplantation. (E and F) Composite overlay of hNCAM+ fibers from DAB-developed sections from 7 animals transplanted with human fetal cell (purple, E) and 9 animals transplanted with VM-patterned hPSC (green, F) intranigral grafts, innervating the host brain. G, H, hNCAM⁺ fibers from DAB-developed sections in G, ventromedial PFC and, (H) dlSTR. (Scale bar, 200 μm.) (I) Immunocytochemistry of TH/ALDH1A1/GIRK2 and (J) TH/CALB in VM-patterned hPSC intranigral grafts. (Scale bar, 50 μm.) (K) Immunohistochemistry of DAT under bright-field illumination of DAB-developed sections. (Scale bar, 50 μm.) (L) Fontana Masson/cresyl violet immunohistochemistry of VM-patterned hPSC intranigral grafts 12 mo posttransplantation. (Scale bar, 20 μm.) (M and N) Uniform manifold approximation and projection (UMAP) embeddings showing clustering of analyzed cells from (M) fetal VM-derived, and (N) VM-patterned hPSC intranigral grafts 12 mo posttransplantation. For both graphs the number of cells were randomly downsampled to 15 000. Cell-type assignments are indicated. (O) Boxplot showing higher hetereogeneity in fetal derived grafts compared to VM-patterned analyzed samples as indicated by a decreased Rogue score. (P) Proportion of each cell type from fetal and VM-patterned hPSC intranigral grafts 12 mo posttransplantation. (Q) UMAP colored by projection of fetal snRNA data onto hPSC-derived cells, colored by prediction score. Blue indicates low similarity (score ~0.4) and brown indicates high similarity (score ~1.0). (R) Violin plot showing the prediction scores across different cell types derived from VM-hPSC. Higher prediction scores, close to 1.0, indicate a strong resemblance to the corresponding fetal cell types, with neurons, OPC, and astrocytes having the highest scores. (S) Dot plot comparing fetal and VM-hPSC-derived cells across various cell types. The size of each dot represents the proportion of cells overlapping.

Fig. 2.

Human DA neuron composition in hPSC-derived grafts. (A) UMAP plot displaying neuron subclusters after reclustering of VM-patterned hPSC intranigral grafts 12 mo posttransplantation (n = 12 509). (B) Violin plots illustrating normalized expression levels of indicated genes in each neuronal subcluster (1 to 7). (C) UMAP plot revealing neuron subclusters after reclustering of fetal-VM derived neurons 12 mo posttransplantation (n = 3954). (D) Heatmap presenting differentially expressed genes and manually selected markers in fetal neuron subclusters (1 to 7). (E) Quantification of neuronal subtype diversity between fetal VM-derived and VM-patterned hPSC intranigral grafts 12 mo posttransplantation (P < 0.01). (F) UMAP of hPSC- derived neuronal cells marked by TH expression (Expression > 0 marked in orange). (G) UMAP plot of DA neuron subclusters after TH⁺ cells reclustering. (H) Dendrogram of unsupervised hierarchical clustering of DA clusters using uncentered Pearson correlation. (I) Individual UMAP plot presenting DA cluster diversity within hPSC-derived intranigral graft from 3 to 12 mo. (J) Heatmap displaying differentially expressed genes and selected markers for DA neuron subclusters (1 to 9). (K) Violin plots showing expression levels of indicated genes in each DA neuron subcluster (1 to 9).

Fig. 3.

Axon target-based single cell transcriptomics of DA neuron subtypes. (A and B) Schematic representation of the TARGET-seq experimental design. Retrograde mCherry AAV injection and tracing from either (A) the ventromedial PFC or (B) the dlSTR to DA neurons in a 12-mo transplant. (C) Detection of retrograde mCherry AAV and TH immunohistochemistry at the ventromedial PFC injection site. (Scale bar, 100 μm.) (D) Detection of retrograde mCherry AAVs, TH immunohistochemistry at the dlSTR injection site. (Scale bar, 100 μm.) (E) Volcano plot illustrating differential gene expression between dlSTR- and PFC-traced human DA neurons. Selected genes are highlighted in green and blue for dlSTR- and PFC-traced human DA neurons, respectively. Dashed lines indicate P-value < 0.05 and |log2-foldchange| > 0.5. (F and G) Gene Ontology (GO) enrichment analysis displaying of the significant genes (P < 0.05) for PFC-traced neurons (F) and dlSTR-traced human (G) DA neurons. Top five enriched ontology terms are shown including leading edge genes for each gene ontology.

Fig. 4.

Gene-set and enrichment pathway analyses in DA subtypes. (A) Cell density plot of PFC- and dlSTR-projected DA neurons in the DA UMAP space. (B) Permutation test of cluster proportions in PFC- and dlSTR-projected DA neurons. Dot-whiskers marked in blue and green are statistically significant (FDR P < 0.05). (C) Violin plots illustrating differential expression levels of indicated genes in dlSTR and PFC -projected DA neuron clusters (2 and 6 respectively). (D) Active regulons defining PFC-projected DA neuron clusters. (E) Gene Set Enrichment Analysis (GSEA) identifying WNT Canonical Signaling in PFC-projected DA cluster (FDR P < 0.05). (F) Heatmap of differentially expressed (P < 0.05) WNT-pathway associated genes between PFC- and dlSTR-projected DA neuron clusters. (G) Active regulons defining dlSTR-projected DA neuron cluster. (H) GSEA identifying MAPK pathway in the dlSTR-projected DA cluster. (I) Heatmap presenting differentially expressed MAPK-pathway associated genes between PFC- and dlSTR-projected DA neuron clusters. (J) MAPK-pathway enrichment analysis in the DA vulnerable population from human postmortem ventral midbrain (VM) compared to all other clusters in the same dataset (13).

Fig. 5.

Novel TFs in Human DA Neuron Subtype Trajectories. (A) Force-directed graph (SPRING) identifying the trajectory of DA cluster in stem cell transplant from month 3 to month 12. Cells are color-graded based on months after grafting. (B) Trajectory reconstruction by minimum spanning tree (Slingshot) of identified DA neuron clusters along the temporal axis. Branches and molecular specification are visualized with a force-directed layout, where each cell is represented by a point proceeding from 3 to 12 mo. (C) Selected markers along the temporal axis showing differential expression along identified trajectories. (D) Gene regulatory network underlying PFC-projected, and (E) dlSTR-projected DA neuron clusters. GRNs were constructed using Celloracle. (F) Quiver plot of CellOracle in silico knockout and overexpression of TFs not previously linked to DAergic neurons with a specific projection pattern: TCF7L and NFIB for PFC-projected DA cluster, and (G) KLF7 and E2F3 for dlSTR-projected DA cluster. The effect of the perturbation is shown with projections of cell state transition vectors on each cell’s UMAP plot. Yellow or red arrow colors indicate in or out directionality, respectively. (H) Representative bright-field image of hPSC-derived DA neurons in 2D culture. (Scale bar ,200 µM.) (I) Cell density plot of projected 2D VM culture DA scRNAseq dataset onto hPSC-derived transplants. (J) Representative bright-field image of VM organoid. (Scale bar, 1 mm.) (K) Cell density plot of projected of VM organoid DA scRNAseq dataset onto hPSC-derived transplants. (L) Representative bright-field image of VM-STR assembloid;VM organoid region within assembloid was labeled by GFP. (Scale bar, 100 µM.) (M) Cell density plot of projected assembloid DA scRNAseq dataset onto hPSC-derived transplants. (N) Quantification of DA developmental trajectories detected in 12 mo hPSC-derived transplants in 2D, VM organoid at different developmental stage (day 60 to 120) datasets. (O) Quantification of PFC- and dlSTR-projected DA cluster in VM organoid and VM-STR assembloid DA datasets.