Single-cell RNA-sequencing of differentiating iPS cells reveals dynamic genetic effects on gene expression

Abstract

Recent developments in stem cell biology have enabled the study of cell fate decisions in early human development that are impossible to study in vivo. However, understanding how development varies across individuals and, in particular, the influence of common genetic variants during this process has not been characterised. Here, we exploit human iPS cell lines from 125 donors, a pooled experimental design, and single-cell RNA-sequencing to study population variation of endoderm differentiation. We identify molecular markers that are predictive of differentiation efficiency of individual lines, and utilise heterogeneity in the genetic background across individuals to map hundreds of expression quantitative trait loci that influence expression dynamically during differentiation and across cellular contexts.

Fig. 1. Single-cell endoderm differentiation of pooled iPSC lines.

a Overview of the experimental design. iPSC lines from 125 genotyped donors were pooled in sets of 4-6, across 28 experiments, followed by differentiation towards definitive endoderm. Cells were sampled every 24 h (Methods) and molecularly profiled using scRNA-seq and FACS. https://github.com/ebiwd/EBI-Icon-fonts by EBI Web Development is licensed under CC BY 4.0. b Variance component analysis of 4,546 highly variable genes, using a linear mixed model fit to individual genes to decompose expression variation into time point of collection, cell line and experimental batch (Methods). c Top: Principal component analysis of gene expression profiles for 36,044 QC-passing cells. Cells are coloured by the time point of collection. Bottom: Cells are ordered by pseudotime, defined as the first principal component (PC1). From left to right, cells transition from a pluripotent state to definitive endoderm. d Single cell expression (y-axis) of selected markers for each developmental stage, spanning iPSC (NANOG), mesendo (T), and defendo (GATA6) stages, plotted along pseudotime (x-axis).

Fig. 2. Mapping single-cell eQTL in each developmental stage.

a Illustration of the single cell eQTL mapping strategy at different stages of differentiation. Shown is an example of an eQTL that is specific to the definitive endoderm (defendo) stage. Boxplots of gene expression stratified by the allelic state of rs9648854 at each stage, showing an association between rs9648854 and CNTNAP2 expression at the defendo stage, but not at earlier stages. https://github.com/ebiwd/EBI-Icon-fonts by EBI Web Development is licensed under CC BY 4.0. b Comparison of eQTL mapping using different strata of all cells. Stage definition based on pseudotime ordering increases the number of detectable eQTL, compared to using the corresponding time point of collection. Bars represent number of eGenes (genes with at least one eQTL, at FDR < 10%). c Proportion of eQTL that are specific to a single stage, shared across two stages, or observed across all stages (sharing defined as a lead eQTL variant at one stage with nominal significant effects P < 0.05 and consistent direction at another stage). d A lead switching event consistent with epigenetic remodelling. The overlap of H3K4me1 with the eQTL SNPs across differentiation time points is indicated by the coloured bars.

Fig. 3. eQTL dynamics during differentiation.

a Combined analysis of the gene expression, ASE, and eQTL dynamics across pseudotime. Upper panels: Schematic of sliding window approach. Cells are binned according to pseudotime groups, to quantify average expression, perform an eQTL analysis, and quantify average ASE (each bin includes 25% of cells, binned at increments of 2.5%). Lower panels: clustered heatmap of expression levels, eQTL effects, and ASE across pseudotime for the top 311 genes with the strongest dynamic QTL effects (FDR < 1%; out of 785 at FDR < 10%; Methods). For each gene, the expression and the ASE dynamics were jointly grouped using clustering analysis, with 4 clusters. The membership of gene expression and ASE dynamics of these 4 clusters is indicated by colours in the right-hand panel. Values in all heatmaps are z-score normalised by row. For ASE, average ASE values are plotted such that red indicates highest deviation from 0.5. b Summary of the identified cluster dynamics, displaying the average dynamic profile of each cluster, computed as the average across z-score normalised gene expression/ASE profiles. c Exemplars of the dynamic gene expression and dynamic genetic effects clusters shown in a. Shaded regions indicate standard error (+/− 1 SEM; Methods). d Number of genes categorised by the combination of expression and ASE cluster from a. Average dynamics of expression clusters (rows) and ASE clusters (columns) as in b are shown. e Overlap of dynamic eQTL variants from a with histone marks. The odds ratio compared to the background of all other eQTL variants is shown (*P < 0.01; **P < 1 × 10⁻⁴; Fisher’s exact test).

Fig. 4. Allele-specific expression reveals interactions with fundamental cellular processes.

a Illustration of eQTL affected by cellular context. Left: Schematic of cellular contexts affecting a regulatory element containing an eQTL SNP, and thus affecting allele-specific expression. Right: Allele-specific expression variation for two exemplar eQTL SNPs that tag cancer GWAS variants and display GxE interactions (FDR < 10%). The eQTL for RNASET2 (rs2247315) tags a risk variant for basal cell carcinoma, and is responsive to cellular respiration, while that for SNRPC (rs9380455) tags a risk variant for prostate cancer and is responsive to the cell cycle G2/M transition (Table S13). Cellular contexts for each cell were inferred by GO annotations of coexpression modules (Methods). Shaded regions indicate standard error (+/− 1 SEM; Methods). b Results summary: numbers of eQTL (from Fig. 2; Methods) identified as displaying GxE interactions with pseudotime (purple), displaying GxE interactions with other cellular contexts but not with pseudotime, (after appropriately accounting for pseudotime, red), displaying GxE interactions with both pseudotime and at least one other cellular context (yellow), and displaying no GxE interactions at all (grey). Significance is assessed at FDR < 10%. c Higher order interaction example: an eQTL variant for EIF5A (rs7503161) is affected by a GxExE higher order interaction with both pseudotime and the G1/S transition. Left panel: Effects of respiration state on ASE for cells with low and high pseudotime. Lines shown are linear regressions with 95% confidence intervals for the 30% of cells with lowest and highest values for pseudotime. Right panel: Heatmap of averaged ASE for cells falling within the specified windows of pseudotime and respiration state. Only values for windows containing n > 30 cells are shown (n = 6423 cells in total).

Fig. 5. Identification of molecular markers for differentiation efficiency.

a Variation in differentiation efficiency across cell lines. Left: Differentiation progress over time, showing trajectories for 98 cell lines, coloured by differentiation efficiency. Shown are 98 cell lines with sufficient data at all time points (out of 126, more than 10 cells). Differentiation efficiency of a cell line was defined as the average pseudotime across all cells on day 3. Right: Differentiation efficiency across cell lines (points), and consistency of individual cell lines differentiated in multiple experiments (vertical bars). b Associations between iPSC gene expression levels and differentiation efficiency. Left: schematic. Centre: Genome-wide analysis to identify markers of differentiation efficiency, considering iPSC gene expression levels. Displayed are negative log P-values signed by the direction of the effect. Horizontal blue lines denote FDR = 10% (Benjamini-Hochberg adjusted). Autosomal genes with significant associations are shown in blue; X chromosome genes with significant associations are shown in red. Right: Scatter plot between gene expression in the iPS state and differentiation efficiency for the X chromosome gene ZDHHC9.