Population-scale single-cell RNA-seq profiling across dopaminergic neuron differentiation

Abstract

Studying the function of common genetic variants in primary human tissues and during development is challenging. To address this, we use an efficient multiplexing strategy to differentiate 215 human induced pluripotent stem cell (iPSC) lines toward a midbrain neural fate, including dopaminergic neurons, and use single-cell RNA sequencing (scRNA-seq) to profile over 1 million cells across three differentiation time points. The proportion of neurons produced by each cell line is highly reproducible and is predictable by robust molecular markers expressed in pluripotent cells. Expression quantitative trait loci (eQTL) were characterized at different stages of neuronal development and in response to rotenone-induced oxidative stress. Of these, 1,284 eQTL colocalize with known neurological trait risk loci, and 46% are not found in the Genotype-Tissue Expression (GTEx) catalog. Our study illustrates how coupling scRNA-seq with long-term iPSC differentiation enables mechanistic studies of human trait-associated genetic variants in otherwise inaccessible cell states.

Ext Fig 1

Ext Fig 2

Ext Fig 3

Ext Fig 4

Ext Fig 5

Ext Fig 6

Ext Fig 7

Ext Fig 8

Fig. 1. Experimental design and cell type heterogeneity in pooled differentiations of iPSCs to a midbrain cell fate.

(a) Experimental workflow for scRNA-seq analysis of iPSC-derived dopaminergic neurons. Time points at which cells were collected for scRNA-seq profiling (day 11, day 30, day 52) are indicated. On day 51, half of the cells were stimulated with rotenone (ROT) for 24 h to induce oxidative stress. (b) UMAP plot of all 1,027,401 cells assayed, colored by annotated cell type identity (Methods). Cells that were not collected at a given condition (time point, stimulus) are displayed in light grey. Prolif: Proliferating. (c) Bar plot showing, for each condition, the fraction of cells assigned to each cell type.

Fig. 2. Reproducible variation in differentiation trajectories.

(a) Box plots display the proportions of each cell type at day 52 across cell lines. The sum of the proportions of DA and Sert cells at day 52 is defined as neuronal differentiation efficiency. In the box plots, the middle line is the median and the lower and upper edges of the box denote the first and third quartiles. (b) Left: hierarchical clustering of (cell line, pool) combinations by neuronal differentiation efficiency. Displayed are data from 10 pools with at least 10 scRNA-seq profiles across all time points (138 lines). First bar indicates a line’s pool of origin; second bar the replicate status. Right: UMAPs, highlighting the distributions of cells on day 52 for two selected cell lines with low (HPSI0514i-fiaj_1; seagreen) and high (HPSI1213i-hehd_1; dark blue) neuronal differentiation efficiencies. (c) Workflow for scRNA-seq analysis of iPSC-derived cerebral organoids. UMAPs summarize the resulting cell populations (1: Floor Plate Progenitors (FPP), 2: Proliferating FPP 3: Neuroblasts, 4: Dopaminergic Neurons, 5: Serotonergic-like neurons (Sert), 6: Proliferating Sert, 7: Astrocyte-like, 8: Ependymal-like (Epen) 1, 9: Epen2, 10: Unknown_Neur1, 11: Unknown_Neur2, 12: Unknown_Neur3; a: Neurons, b: Intermediate progenitors, c: Radial glial progenitors, d: Satellite cells, e: Mesenchymal cells, f: Myotube, g: PAX7⁺ cells, h: Wnt⁺ cells). (d) Scatter plot of neuronal differentiation efficiency between replicate pools (n = 32 cell lines differentiated in two different pools). Highlighted are the two cell lines from b. (e) UMAPs of the two cell lines selected in b, making non-brain and brain cell types in the organoid study. (f) Scatter plot of midbrain dopaminergic neuronal differentiation (x-axis) versus neuronal differentiation efficiency measured in organoid differentiation (y-axis), for a common subset of 12 iPSC lines. Highlighted are the two cell lines from b. In panels d and f, LOESS curves and 95% confidence intervals alongside Pearson’s R and the P value from a two-sided t-test.

Fig. 3. A gene expression signature in iPSCs is associated with neuronal differentiation efficiency.

(a) Venn diagram indicating the overlap of cell lines included in this study and two recent iPSC studies, a bulk RNA-seq study and a single-cell RNA-seq study. (b) Histogram of Pearson correlation coefficients between variation in gene expression of individual genes (from bulk RNA-seq) and neuronal differentiation efficiency. Two representative genes with positive (TAC3) and negative (UTF1) association are highlighted. Scatter plot of neuronal differentiation efficiency (x-axis) versus bulk iPSC gene expression (y-axis) for UTF1 and TAC3. (d) UMAPs of single-cell RNA-seq profiles in iPSCs from 112 donors. Colors denote the expression level of the two selected genes from b,c: UTF1 and TAC3. Cluster 2 is shown by the dashed lines. (e) Comparison of marker gene association results with expression markers of the cluster 2. Shown is, for each gene, a scatterplot of the Pearson correlation coefficient of association between iPSC gene expression and neuronal differentiation efficiency (x-axis, based on bulk gene expression) versus its log fold-change between cluster 2 and all other clusters (y-axis, based on scRNA-seq). The genes UTF1 and TAC3 are highlighted. (f) Scatter plot between neuronal differentiation efficiency (x-axis) and the proportion of cells assigned to cluster 2 (y-axis) across 45 cell lines that were included in both sets of experiments. Where multiple measurements were available for a given cell line, average values are shown. In panels c, e and f, the Pearson’s R and the P value from a two-sided t-test are indicated. Diff. = differentiation.

Fig. 4. Mapping of cis eQTL in 14 distinct cell contexts (“cell types”-”conditions”) for 6 dominant cell types identified across midbrain differentiation.

(a) Cumulative number of genes with at least one eQTL (eGenes) for each cell type and condition (D11 = day 11; D30 = day 30; D52 = day 52; ROT = Rotenone stimulation). (b) Left: day 52-specific eQTL for HSPB1 in DA (rs6465098; FDR <5%, Methods). Shown are Manhattan plots for DA at day 30 (top) and unstimulated DA at day 52 (bottom). Right: a rotenone stimulus-specific eQTL for ACSF3 in Sert (rs12597281, right). Shown are Manhattan plots for rotenone-stimulated Sert at day 52 (top) and unstimulated day 52 Sert (bottom). (c) Comparison of the number of eGenes for individual eQTL maps (FDR <5%; y-axis) as a function of the effective sample size (number of unique donors; x-axis) across studies and cell types. Left: results from overlapping eQTL results in this study with in vivo eQTL maps from GTEx, divided into brain tissues (yellow) and non-brain tissues (grey). In red, the result from this study when aggregating across all 14 eQTL maps. Right panel shows a magnified view of results from our study colored by cell type and shaped by condition. (d) Sharing of eQTL signals discovered across 14 cell contexts in this study, as well as using bulk RNA-seq in iPSCs^,, with in vivo brain eQTL maps (from GTEx). Violin plots show the extent of eQTL sharing (Methods), with each of 13 GTEx brain eQTL maps. In the box plots, the middle line is the median and the lower and upper edges of the box denote the first and third quartiles, while the violin plots show the distribution. Astro: Astrocyte-like; DA: Dopaminergic neurons, Epen1: Ependymal-like1, FPP: Floor Plate Progenitors, P_FPP: Proliferating FPP, Sert: Serotonergic-like neurons, U_Neur1,2,3: Unknown neurons 1,2,3.

Fig. 5. Colocalization analysis of eQTL with 25 neuro-related GWAS traits.

(a) Venn diagram showing the number and the overlap of colocalization events discovered using eQTL maps from this study, GTEx brain and other GTEx tissues. (b) Heatmap showing the posterior probability of colocalization (PP4 from COLOC; Methods) for eQTL that colocalized with one or more GWAS traits. N: Neuronal differentiation (this study), B: GTEx Brain, O: Other GTEx tissues. (c) Locus zoom plots around the SFXN5 gene. The schizophrenia GWAS association (left) is colocalized with the eQTL in rotenone-stimulated serotonergic-like neurons at D52 (second panel from the left). No colocalization signal was detected in unstimulated serotonergic neurons at D52 (third panel from the left) or any other brain GTEx tissues as illustrated here with GTEx brain amygdala (rightmost panel). The lead variant is indicated with a purple diamond and other points were colored according to the LD index (r2 value) with the lead variant. (d) A midbrain progenitor-specific eQTL for FGR1 associated with schizophrenia. We identified a colocalization event with this eQTL in both proliferating (second panel from the left) and non-proliferating floor plate progenitors (third panel from the left) at day 11. No colocalization was found in any other cell type profiled in this study (not shown) nor in any brain GTEx tissues (shown with GTEx brain hypothalamus, rightmost panel). Astro: Astrocyte-like; DA: Dopaminergic neurons, Epen1: Ependymal-like1, FPP: Floor Plate Progenitors, P_FPP: Proliferating FPP, Sert: Serotonergic-like neurons, U_Neur1,2,3: Unknown neurons 1,2,3.