Molecular and functional variation in iPSC-derived sensory neurons

Abstract

Induced pluripotent stem cells (iPSCs), and cells derived from them, have become key tools for modeling biological processes, particularly in cell types that are difficult to obtain from living donors. Here we present a map of regulatory variants in iPSC-derived neurons, based on 123 differentiations of iPSCs to a sensory neuronal fate. Gene expression was more variable across cultures than in primary dorsal root ganglion, particularly for genes related to nervous system development. Using single-cell RNA-sequencing, we found that the number of neuronal versus contaminating cells was influenced by iPSC culture conditions before differentiation. Despite high differentiation-induced variability, our allele-specific method detected thousands of quantitative trait loci (QTLs) that influenced gene expression, chromatin accessibility, and RNA splicing. On the basis of these detected QTLs, we estimate that recall-by-genotype studies that use iPSC-derived cells will require cells from at least 20-80 individuals to detect the effects of regulatory variants with moderately large effect sizes.

Figure 1. Characterization of molecular phenotypes in iPSC-derived sensory neurons.

(a) Schematic of IPSDSN differentiation and assays. iPSCs were received in Essential 8 (E8) medium (N=82) or on mouse embryonic fibroblasts (MEFs, N=49), and transferred to KSR-XF medium. Over 11 days, different inhibitor combinations were added (2i, 5i, 3i, see Methods), and N2B27 medium phased in, followed by transfer to growth factor medium at day 11 for neuronal maturation. (b) PCA plot projecting IPSDSN, iPSC, and DRG samples onto the first two principal components defined based on RNA-seq FPKMs in GTEx tissues. Some GTEx tissues are unlabelled due to overlapping labels. (c) Expression of sensory neuronal marker genes (SCN9A, DRGX) and key iPSC genes (NANOG, POU5F1).

Figure 2. Single-cell sequencing of IPSDSN cells.

(a) A heatmap of RNA-seq data for ten marker genes of the two cell clusters identified by SC3. Color scale denotes normalised gene expression levels. (b) The first two principal components (PCs) of IPSDSN gene expression, with estimated fibroblast-like percentage from CIBERSORT, from samples derived using protocols 1 and 2 (P1 and P2).

Figure 3. Gene expression variability in IPSDSNs is influenced by differentiation conditions.

(a) Density plot of the coefficient of variation of genes across samples, separately for each GTEx tissue, IPSDSN samples (n=106, P2 protocol only), iPSC (n=239), and DRG (n=28). (b) Violin plot showing, for each gene, the estimated fraction of total expression variability across samples due to differentiation batch, donor genetics or iPSC reprogramming, culture conditions (“wasFeeder”: feeder-dependent vs. E8 medium), and gender. (c) Differentially expressed genes (FDR 1%, blue and red points) between iPSC samples grown on feeders (n=68) vs. E8 medium (n=171). (d) Differentially expressed genes (FDR 1%) between IPSDSNs from feeder- (n=27) and E8-iPSCs (n=79). Neuronal differentiation genes, such as RET and L1CAM, are more highly expressed in samples from E8-iPSCs. (e) Left boxplot: global gene expression differences between feeder- and E8-iPSCs are captured in PC1. Right two boxplots: selected differentially expressed genes. (f) Left boxplot: estimated neural fraction of samples differs in IPSDSNs derived from feeder- and E8-iPSCs. Right two boxplots: selected differentially expressed genes. Boxplots show the median, 25th and 75th percentiles, with whiskers extending 1.5 times the interquartile range.

Figure 4. Splicing QTLs overlapping GWAS.

(a) An sQTL for TNFRSF1A leads to skipping of exon 6, and overlaps with a multiple sclerosis association. (b) An sQTL for SIPA1L2 leads to increased skipping of an unannotated exon between alternative promoters, and overlaps with a Parkinson’s disease association. (c) An sQTL for APOPT1 alters skipping of exons 2 and 3, and overlaps with a schizophrenia association. P values are from the FastQTL beta approximation based on 10,000 permutations. Boxplots show the median, 25th and 75th percentiles, with whiskers extending 1.5 times the interquartile range.

Figure 5. Power to detect a genetic effect in a single-variant single-gene test depends on sample size, allelic effect size, and gene expression variability.

(a) TPR as a function of allelic fold change for five different numbers of replicates (half the total sample size). (b) TPR as a function of CV for five bins of allelic fold change, with 10 samples of each genotype.