Differentiation progress and cell line differences drive variation in gene expression

In order to characterize the complex landscape of cardiomyocyte differentiation, we used UMAP to produce a low-dimensional embedding of the single cell data while preserving global structure. We found that while cells from the early days of the differentiation time course exhibited fairly uniform transcription profiles, this was less true for later days (days 7, 11, and 15; Fig 1A and 1D). Marker genes known to be expressed at various stages in cardiac differentiation, from iPSC to mesoderm to cardiomyocyte, showed high expression at expected early, intermediate, and late stages of the differentiation time course, respectively (Fig 1B and 1E). Next, we used unsupervised clustering to partition cells into clusters, and matched cell clusters to known cell types based on expression of known marker genes (Fig 1C, Materials and Methods, [36]). While these discrete cell type assignments are imperfect as they fail to capture the continuity of the differentiation process, they are useful in characterizing the broad relationships between groups of cells. As suggested by previous reports [23,35], we identified a bifurcation in the differentiation landscape, giving rise to two distinct terminal cell types. One of these terminal cell types has high expression of genes known to be involved in cardiomyocyte function, such as TNNT2 and MYL7 (Fig 1B, [37,38]). Cells in the other terminal cell type do not express cardiomyocyte markers, and instead have high expression of genes such as COL3A1 and VIM, which are expressed in the extracellular matrix of cardiac fibroblasts [39,40]. The differentiation outcome of each sample, namely the proportion of cells in each cluster, varied by individual cell line; certain lines differentiated primarily into either the TNNT2-expressing or the COL3A1-expressing terminal cell type clusters (Figs 1E and S3). For the remainder of this paper, we will refer to the TNNT2-expressing cell cluster as cardiomyocytes (CM) and to the COL3A1-expressing cluster as cardiac-fibroblasts (CF) or fibroblast. We also identified a cluster that underexpressed marker genes of cardiac cell types throughout the differentiation process, and instead expressed several endoderm-specific markers such as APOA1 and AFP. We were unable to fully characterize this cluster based on expression patterns alone, and omitted these cells from downstream investigation of the dynamics of gene regulation on gene expression during mesoderm and cardiac cellular differentiation (S4 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Gene expression patterns in single cell data.

(A) UMAP of full single cell dataset; cells are colored by differentiation day. (B) Estimated density of expression for several marker genes across cells. (C) UMAP of full single cell dataset; cells are colored by cell type, assigned based on Leiden clustering and marker gene expression. IPSC = induced pluripotent stem cell, MES = mesoderm, CMES = cardiac mesoderm, PROG = cardiac progenitor, CM = cardiomyocyte, CF = cardiac fibroblast, UNK = unknown cell type. (D) Proportion of cells belonging to each cell type per differentiation day, across all cell lines. (E) Distribution of L1TD1 (pluripotency marker), TNNT2 (cardiomyocyte marker) and COL3A1 (cardiac fibroblast marker) over cells from 6 representative examples of the 19 cell lines studied, for each of the 7 differentiation days.

https://doi.org/10.1371/journal.pgen.1009666.g001

Single-cell expression data resolves bifurcating trajectories during cellular differentiation

In previous work, we investigated the relationship between genotype and chronological time, represented by the differentiation day in which each bulk sample was collected. However, chronological time may not properly capture the axis of variation along which genetic regulation is changing, and can be heavily confounded by heterogeneity in differentiation within and between samples. If cells within a sample progress through differentiation at different rates, their aggregated expression profile will not be truly reflective of an individual stage of differentiation, confounding tests for association between genotype and differentiation progress. Systematic differences between cell lines can exaggerate this: differentiation speed appears to vary between cell lines (Fig 1E), such that differentiation progress at day 3, for example, is not uniform across samples. Such differences can lead to false associations between genotype and differentiation progress in cases where genotype is partially correlated with a cell line’s differentiation speed.

Cellular heterogeneity drives further challenges when aggregating across cells that are differentiating along diverging paths. Aggregated bulk profiles will lose information about the individual cell types present, and if cell type composition varies between individuals (Fig 1E), this will further confound associations between genotype and expression changes during differentiation.

By collecting expression at the single-cell level, we are able to address both of these challenges. To properly focus on the two primary cardiac lineages present, we used the scanpy package to produce a low-dimensional Force Atlas embedding of the cells that had been successfully assigned to a known cell type (Fig 2A, Materials and Methods, [41,42]). We inferred pseudotime for each cell with diffusion pseudotime [43,44], so that progress through differentiation is learned from cells’ individual expression profiles rather than their time of collection (Fig 2B). We performed trajectory analysis using PAGA [44] to examine the relationships between the cell types present in this dataset, which helped to resolve two distinct lineages present in the data, giving rise to cardiomyocyte and cardiac fibroblast cell types (S5 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Pseudotime inference and pseudobulk aggregation.

(A) Force atlas embedding of all cells from the two cardiac differentiation lineages combined, colored by cell type. IPSC = induced pluripotent stem cell, MES = mesoderm, CMES = cardiac mesoderm, PROG = cardiac progenitor, CM = cardiomyocyte, CF = cardiac fibroblast, UNK = unknown cell type. (B) Force atlas embedding from (A), colored by log₂(pseudotime+1), which was inferred for each cell shown using diffusion pseudotime. (C) Distribution of normalized expression of TNNT2, a cardiomyocyte marker gene, across cells from the cardiomyocyte lineage for each differentiation day. (D) Normalized TNNT2 expression across cells from each of 16 pseudotime quantile bins along the cardiomyocyte trajectory. (E) Normalized expression of COL3A1, a cardiac fibroblast marker, across cells from the cardiac fibroblast lineage for each differentiation day. (F) COL3A1 expression across cells for 16 pseudotime quantile bins along the cardiac fibroblast trajectory.

https://doi.org/10.1371/journal.pgen.1009666.g002

One disadvantage to single-cell data compared to bulk is that single-cell measurements are more sparse and noisy: by aggregating over cells, bulk RNA-sequencing reduces noise, which makes expression measurements more tractable for eQTL calling. This introduces a tradeoff between the flexibility of analysis at the single cell level, where we can explore a broader range of dynamic effects among finely resolved pseudotimes or cell populations, and the robustness of analysis on aggregated data [45]. To balance this tradeoff, we partitioned cells from each lineage into pseudotime bins, pooling information across cells to mitigate the noisiness of single cell expression measurement while maintaining homogeneous populations of cells through lineage subsetting and pseudotime binning. Cells were assigned based on the trajectory analysis to the cardiomyocyte lineage, the cardiac fibroblast lineage, or both in the case of precursor cell types. This aggregation scheme enables us to produce a greater number of samples, as we are no longer constrained to the 7 days when experimental collection was performed, while maintaining the expected trends of lineage-specific marker gene expression over pseudotime (Fig 2C-2F).

Mapping of dynamic eQTLs

We applied a Gaussian linear model to the aggregated single-cell pseudo-bulk data based on pseudotime bins from each lineage to identify dynamic eQTLs, namely variant-gene pairs in which the interaction effect of genotype and differentiation time is significantly associated with changes in gene expression. We identified linear dynamic eQTLs for 357 genes in the cardiomyocyte lineage (q < 0.05) and 903 genes in the cardiac fibroblast lineage (Materials and Methods; Table 1). The difference in the number of dynamic eQTLs detected between the two lineages may arise due to greater heterogeneity among predominantly cardiomyocyte samples, or from a difference in the number of cells captured (45,980 cardiac fibroblasts versus 21,862 cardiomyocytes) leading to more precise pseudobulk expression profiles in predominantly cardiac fibroblast samples. We found that both lineage specificity and the replacement of real chronological time with pseudotime improved power for dynamic eQTL detection: using chronological differentiation day as the time variable identified only 142 and 29 dynamic eQTLs for the cardiomyocyte and cardiac fibroblast lineages, respectively. Using differentiation day as the time variable and omitting lineage specificity altogether identified only 5 dynamic eQTLs in the pseudobulk data. For comparison, an analysis of static (non-dynamic) eQTLs using the mashr framework [46], which similarly aggregated pseudobulk by cell line and differentiation day, revealed 183 static eQTLs (Materials and Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Comparison of linear dynamic eQTL calling methods.

We report the number of linear dynamic eGenes (genes with a significant dynamic eQTL at gene-level q-value < = 0.05), for each of the aggregation schemes assessed. Total number of genes tested and total number of tests run are also reported.

https://doi.org/10.1371/journal.pgen.1009666.t001

Ultimately, our lineage subsetting and pseudotime approach revealed slightly more dynamic eQTLs than were previously identified in an experiment with bulk collections at over twice as many time points [23]. To ensure a meaningful comparison, we re-processed the previously collected bulk data in a similar pipeline as pseudo-bulk, accounting for changes in hypothesis testing and filtering of variant-gene pairs (Materials and Methods). This revealed a total of 1028 genes with a dynamic eQTL (compared to a total of 1056 genes detected between both lineages with pseudobulk binned to a similar number of samples; S6 Fig). When the bulk data was subset to the same 7 collection time points used for the single cell experiment, only 210 dynamic eGenes were detected. The increased detection rate offered by the pseudobulk analysis may stem from increased homogeneity of cellular populations that undergo pseudobulk aggregation, as well as improved measurement of differentiation progress achieved by using cellular pseudotime rather than sample collection time.

We investigated several potential sources of confounding in this analysis. First, we assessed calibration of dynamic eQTL calling by permuting the pseudotime variable in the interaction term before calling dynamic eQTLs, which did not reveal evidence of substantial inflation (Materials and Methods, S7 Fig). Second, we investigated the possibility that broad differences between cell lines, such as variation in differentiation speed or trajectory preference, are driving false positive discoveries. We did not find evidence of confounding among dynamic eQTLs, as would be suggested by elevated pairwise correlation of genotype at these loci (Materials and Methods, S8 Fig). This is likely because regression of cell line principal components (Materials and Methods, [23]) effectively controls for such broad differences (S9 Fig). Third, we checked for type I error inflation due to ’double dipping’, where the hypotheses tested are influenced by analysis of the data itself [47–49], since pseudotime is inferred from single-cell expression data. We demonstrate in simulation that the fixed-effect linear model used in this study was conservative in the presence of multiple measurements per individual and did not lead to type I error inflation, though there may be some loss in sensitivity (Materials and Methods, S10 Fig). Finally, we investigated the impact of uncertainty in pseudotime estimation by examining the impact of permuting cellular pseudotime before pseudobulk aggregation, which did not suggest inflation (Materials and Methods, S11 Fig).

As an example of the trait relevance of these dynamic eQTLs, one dynamic eQTL variant, rs1234988, has previously been implicated by GWAS to be associated with hypertension (p = 2.5e-35), and was detected as a dynamic eQTL for ARHGAP42, a Rho GTPase which has previously been identified as a critical regulator of vascular tone and hypertension in mice (Figs 3A, 3B and S12, [50,51]). Notably, ARHGAP42 is known to be a smooth-muscle selective Rho GAP, and this dynamic eQTL was exclusively identified in the cardiac fibroblast lineage (Bonferroni-adj. p = 2.4e-5, cardiac fibroblast lineage, adj. p = 0.79, cardiomyocyte lineage). This variant is not detected as a dynamic eQTL without lineage subsetting or pseudotime binning (adj. p = 1). This example illustrates the advantages of incorporating exploratory data analysis in the study of in vitro experimental datasets: while the differentiation procedure used for these experiments was designed to produce exclusively cardiomyocytes, an alternative terminal cell type discovered after exploratory data analysis is able to provide meaningful insight into an additional differentiation process.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Linear and nonlinear dynamic eQTLs.

(A) rs1234988 is a linear dynamic eQTL for ARHGAP42; the effect of genotype (color) on ARHGAP42 expression (y-axis) varies across pseudotime (x-axis). (B) A previously reported genome-wide association study (bottom) showed that hypertension is associated with genotype at the rs1234988 locus, where a dynamic eQTL for ARHGAP42 was identified. (C) rs1814432 is a nonlinear dynamic eQTL for the gene CFC1B.

https://doi.org/10.1371/journal.pgen.1009666.g003

The pseudotime values can be interpreted as intermediate time points with greater resolution than chronological time. We therefore used these values to also identify nonlinear dynamic eQTLs, whose effects vary in a nonlinear way over the course of differentiation, such as presence only at intermediate stages of the differentiation (Figs 3C and S12). We identified 74 nonlinear dynamic eQTL variants for the cardiomyocyte lineage (q<0.05), and 147 for the cardiac fibroblast lineage (Materials and Methods; S2 Table). Our time course study design is particularly useful for detecting transient nonlinear genetic effects which may not be found by studying only the initial or terminal cell types of a dynamic process such as differentiation.

We examined the extent to which the dynamic eQTLs detected overlapped with eQTLs previously identified in GTEx [4]. After subsetting to gene-variant pairs that were tested in both our data and GTEx, we found that the greatest replication of pseudotime-binned, cardiomyocyte lineage linear dynamic eQTLs occurred in atrial appendage tissue (π₁ = 0.50, method described in [52]), while the greatest replication of pseudotime-binned, cardiac fibroblast linear dynamic eQTLs (as well as bulk) occurred in cultured fibroblasts (π₁ = 0.47, 0.56 respectively). However, by searching directly for dynamic effects across cell types rather than within a single tissue in isolation, we additionally identify eQTLs which were not found to be a significant eQTL in any tissue in GTEx. After subsetting to variant-gene pairs that were tested in both our data and GTEx, we found that 100 of the 359 (28%) linear dynamic eQTLs in the cardiomyocyte lineage were identified as eQTLs in GTEx. Similarly, only 22 of 75 (29.3%) nonlinear dynamic eQTLs on the cardiomyocyte lineage where previously identified as eQTLs in GTEx. Further classification of these cardiomyocyte lineage nonlinear dynamic eQTLs as "early-acting", "late-acting", "transient" (or "middle"), and "switch" (as in ref. [23]) found that most of these variants (68 of 74) are late-acting, supporting the similar replication rates between linear and nonlinear dynamic eQTLs. Of the five early-acting and transient dynamic eQTLs, none were identified as cis-eQTLs in GTEx.

Deconvolution of bulk RNA sequencing data assigns lineage specificity to dynamic eQTLs

The information about the landscape of cardiomyocyte information obtained through single-cell RNA sequencing can also be applied to improve dynamic eQTL calling in bulk data. For each cell type that we identified in the single cell data, we computed a signature expression profile across the top 300 differentially expressed genes that were also measured in bulk (Materials and Methods). We then used CIBERSORTx to deconvolve our bulk data, assigning to each bulk sample a vector of cell type proportions (Figs 4A and S13, Materials and Methods, [53]). Deconvolution reveals that cell type heterogeneity is prominent between samples, particularly in days 7–15. This heterogeneity emphasizes the need to account for cell type proportion in measuring genetic regulatory dynamics, as these broad differences between cell lines can drive false positive associations between time and any genotype that is correlated with broad cell type proportion differences between cell lines.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Cell type deconvolution and interaction eQTL calling.

(A) Cell type deconvolution was applied to decompose RNA expression of a mixed sample, aggregated over multiple cell types, into its constituent cell type proportions (Materials and Methods). Each row represents a cell line, collected in two separate experiments. In the left column, bulk RNA-sequencing data was collected for 15 timepoints (time on x-axis). In the right column, pseudobulk was aggregated across cells collected for 7 time points (time on x-axis). For pseudobulk data, deconvolution is not needed, as each cell is assigned to a cell type. Thus, "ground truth" cell type fractions are accessible as reflected here. (B) Number of genes with a cell type interaction eQTL in bulk for each of six cell types. (C-D) CMYA5 has an interaction eQTL for the cardiomyocyte lineage (C) that is not identified in the cardiac fibroblast lineage (D).

https://doi.org/10.1371/journal.pgen.1009666.g004

We then used these cell type proportions to identify cell type specific regulatory effects, based on cell type interaction eQTLs (ieQTLs) for each known cell type that was observed in the single cell data (Fig 4B). Interaction eQTLs for a cell type at an endpoint of the differentiation (iPSC, cardiomyocyte [CM], and cardiac fibroblast [CF]) are analogous to linear dynamic eQTLs, since cellular composition often partially reflects differentiation time. However, this relationship varies between lineages: we find that cardiomyocyte proportion is more correlated with differentiation day than cardiac fibroblast proportion (Pearson’s ⍴ = 0.59 for cardiomyocyte proportion, compared to ⍴ = 0.36 for cardiac fibroblast proportion). Accordingly, we found that CM and CF ieQTLs called with this approach were replicated in the previously used dynamic eQTL calling framework on the same bulk dataset to varying degrees (π₁ = 0.84 and 0.43, respectively). This may reflect a fundamental difference between the two differentiation trajectories, if cells within samples which produce primarily cardiac fibroblasts reach maturity more quickly. This hypothesis is supported by the observation that the maximum pseudotime value of the cardiomyocyte lineage is greater than that for the cardiac fibroblast lineage (Fig 2B).

This approach enables insights that could not be attained using bulk data alone. The use of an expression signature matrix derived from single cell data with multiple terminal cell types allows us to characterize hundreds of lineage-specific effects which were previously obscured in bulk data: we find that many of the CM- and CF-ieQTLs are lineage-specific (78% and 77% respectively), including some which are potentially relevant to heart-related disease. Fig 4C and 4D show an example of a cardiomyocyte interaction eQTL for cardiomyopathy-associated protein 5 (CMYA5), a gene which is highly expressed in heart and skeletal muscle and has previously been associated with cardiac hypertrophy [54]. This variant was not previously identified by GTEx as an eQTL for CMYA5. More broadly, these interaction eQTLs showed enrichment for genes related to myogenesis that had not been observed among bulk dynamic eQTLs (p = 7e-4, both CM and CF ieQTL, compared to p = 0.17, bulk dynamic eQTL).

This paradigm for analysis is potentially powerful, for example, given bulk data with larger sample sizes or denser time point sampling possible due to lower cost and effort of bulk sequencing, combined with smaller-scale single cell data that offers cell type, lineage, and pseudotime resolution. Given the denser time-point sampling in our bulk data, there are many interaction QTLs that were not detected in the single-cell data alone, or bulk data alone, but enabled by the integration of bulk and single-cell data together.

Ethics statement

The cell lines used in this study were obtained from the NHGRI Sample Repository for Human Genetic Research at the Coriell Institute for Medical Research. All samples were collected by the Coriell Institute for Medical Research with written informed consent and with IRB approval. The genetic data used here has previously been made available through the International HapMap Project [57].

Samples

We used induced pluripotent stem cell (iPSC) lines from 19 individuals from the Yoruba HapMap population. These iPSC lines were reprogrammed from lymphoblastoid cell lines and characterized previously [58]. All 19 individuals were female and unrelated. We chose to use only female individuals to avoid introducing additional variance that is not of interest in this study.

iPSC maintenance

Feeder-free iPSC cultures were maintained on Matrigel Growth Factor Reduced Matrix (CB40230, Thermo Fisher Scientific, Waltham, MA) with Essential 8 Medium (A1517001, Thermo Fisher Scientific) and Penicillin/Streptomycin (30002Cl, Corning, Corning, NY). Cells were grown in an incubator at 37°C, 5% CO2, and atmospheric O2. Cells were passaged to a new dish every 3–5 days using a dissociation reagent (0.5 mM EDTA, 300 mM NaCl in PBS) and seeded with ROCK inhibitor Y-27632 (ab120129, Abcam, Cambridge, UK).

Cardiomyocyte differentiation

We differentiated iPSCs using a protocol previously optimized for use with the Yoruba HapMap panel [58]. This protocol implements slight modifications to the cardiomyocyte differentiation protocols from [59] and [36]. Feeder-free iPSCs were seeded onto wells of a 6-well plate and grown for 3–5 days prior to differentiation. When most lines were 70%-100% confluent, E8 media was replaced with “heart media” along with 1:100 Matrigel hESC-qualified Matrix (08-774-552, Corning) and 12uM of GSK-3 inhibitor CHIR99021 trihydrochloride (4953, Tocris, Bristol, UK). “Heart media” is composed of RPMI (15-040-CM, Thermo Fisher Scientific) with B27 Supplement minus insulin (A1895601, Thermo Fisher Scientific), 2mM GlutaMAX (35050–061, Thermo Fisher Scientific), and 100mg/mL Penicillin/Streptomycin (30002Cl, Corning). CHIR99021 is a small molecule that activates WNT signaling and initiates the differentiation on day 0 (after the ‘day 0’ cell collection) [59]. “Heart media” was replaced 24 hours later at day 1 of differentiation. 48 hours later, at day 3 of differentiation, cells were fed with new “heart media” containing 2uM of the WNT inhibitor Wnt-C59 (5148, Tocris) [59]. We cultured cells in Wnt-C59 heart media for 48 hours. At day 5, Wnt-C59 was removed, and base “heart media” was added. “Heart media” was refreshed on days 7, 10, 12, and 14 of differentiation. Cells began spontaneous mechanical beating between days 7 and 13 of differentiation.

In some cases, after performing cardiac differentiation, one might choose to perform a post hoc purification process to remove any non-cardiac cell types present at the terminal time point [60]. However, for the purposes of a time course experiment where multiple intermediate time points are assayed, a purification protocol undertaken only at the end of the differentiation would not prove useful; therefore, no cell type purification was performed.

Sample collection and processing

We performed cardiomyocyte differentiations in three total batches of six to seven cell lines at a time. For each batch, cardiomyocyte differentiations were performed with three staggered starting days, such that samples could be collected from each cell line in three differentiation stages at any given time. For all 19 cell lines, samples were collected on differentiation days 0 (iPSC, before treatment with CHIR99021), 1, 3, 5, 7, 11, and 15. Drop-seq collection was performed a total of three collection days for each batch of six to seven cell lines. In the first collection day, samples from all cell lines in the batch were collected for differentiation days 1, 3, and 7. In the second collection day, samples from all cell lines in the batch were collected for differentiation days 5, 7, and 11. In the third collection day, samples from all cell lines in the batch were collected for differentiation days 0 (iPSC), 11, and 15. Through this process, single-cell gene expression data was collected for all cell lines in seven unique time points, with two time points (differentiation days 7 and 11) having two replicates. This staggered differentiation and collection study design was performed to minimize the technical effect of sample collection as a potential confounding variable associated with cell line or differentiation day.

To harvest the samples at the start of each collection day, cells in at least two wells of a 6-well culture dish were released from the dish using Accutase (BD Biosciences, San Jose, CA, #561527). Samples were washed three times and resuspended in 1X PBS, 0.01% BSA. Cells were then passed through a 40 um filter to encourage the formation of a single cell suspension. The concentration of each single cell suspension was quantified manually using an NI hemocytometer (INCYTO, Cheonan, Korea, DHC-N01-2).

Using a 125 um Drop-seq microfluidic device, single cells were captured in droplets along with a DNA barcoded bead (ChemGenes, Wilmington, MA, Macosko-2011-10(V+)), following the standard Drop-seq protocol [61]. The DNA barcoded beads include a cell-specific barcode so the cell identity of each RNA molecule can be recovered. After Drop-seq collection, the RNA molecules were reverse transcribed, and cDNA amplification was performed according to the Drop-seq protocol. cDNA concentration and library size were measured using the Qubit 3 fluorometer (Thermo Fisher) and BioAnalyzer High Sensitivity Chip (Agilent, Santa Clara, CA, #5067–4626).

Library preparation was performed using the Illumina Nextera XT DNA Library Preparation Kit (Illumina, FC-131-1096). Libraries in each batch were multiplexed together so that every sequencing lane contained three samples, one from each of the three collection days. Each of those samples was itself a multiplexed collection of three individual cell lines at three distinct differentiation time points, which were mixed upon Drop-seq collection. Samples went through paired-end sequencing using the Illumina NextSeq 500. 20 bp were sequenced for Read 1, and 60 bp for Read 2 using Custom Read 1 primer, GCCTGTCCGCGGAAGCAGTGGTATCAACGCAGAGTAC, according to manufacturer’s instructions [61]. The same multiplexed library pool was sequenced twice with the goal of achieving at least 20 million reads per sample.

We recorded 20 technical and biological covariates and measured their contribution to variation in our data (S14 Fig).

RNA-seq quantification

For each sequencing run, we obtained paired-end reads, with one pair representing the cell-specific barcode and unique molecular identifier (UMI), and the second pair representing a 60 bp mRNA fragment. We used dropseqRunner (available at github.com/aselewa/dropseqRunner) which takes a fastq file with paired-end reads as input and produces an expression matrix corresponding to the UMI of each gene in each cell. All RNA-seq samples were aligned to the human genome (GRCh38) using STAR-solo [62]. We used featureCounts [63] to assign each aligned read to a genomic feature, and umi_tools [64] to create a count matrix representing the frequency of each feature in our dataset. We then used the single-cell demultiplexing software ’demuxlet’ to assign to each cell a probability that the cell is a doublet [65].

The following filter was applied to remove 21,725 rare genes (out of 60,668) from downstream analysis:

• Gene must be detected in at least 10 cells

The following filters were then applied to remove 330,750 low-quality cells (out of 564,362) for downstream analysis:

• Maximum doublet probability of 0.3 from demuxlet
• Unambiguous assignment of the cell to an individual by demuxlet (maintain cells not assigned to ’doublet_ambiguous’)
• Maximum of 25% mitochondrial reads
• Minimum of 300 unique genes detected (of the genes that passed the previous filtering step)

Following these filtering steps, an additional 2,826 cells were removed whose feature or read counts were more than 4 standard deviations away from the median. This left a total of 230,786 cells and 38,943 genes for downstream analysis (S15–S17 Figs).

Cell cycle correction and normalization of single-cell expression data with Seurat

We used the Seurat workflow for cell cycle regression in differentiating. Each cell was assigned a score for G2/M phase and S phase according to marker gene expression, and the difference between these scores was regressed out during normalization. The data was then normalized using the SCTransform function in [66,67], producing corrected counts, log-normalized corrected counts, Pearson residuals, and a set of highly variable features. The Pearson residuals of 1,000 highly variable features were scaled so that each gene had unit variance across all cells for downstream analysis.

Dimensionality reduction and clustering with scanpy

Dimensionality reduction, clustering and pseudotime were performed using the scanpy package [41], following Seurat object to h5ad conversion via the sceasy package [68]. The scaled Pearson residuals from 1000 highly variable features were used to compute 50 principal components (PCs), which were then embedded into a 2D UMAP plot (Fig 1A and 1C). These 50 PCs were also used to produce a neighborhood graph, and Leiden clustering was performed at resolution 0.35 to produce the clusters shown in Fig 1C. (Several clusters are merged into the unknown cell type, as described below).

Lineage specification and pseudotime inference

Based on marker gene expression patterns (Fig 1B), 6 of the 10 Leiden clusters were annotated with known cell types. To facilitate trajectory reconstruction, 3 outlier clusters with less than 5,000 cells were removed. Cluster 7 contained a group of cells which did not express marker genes for cardiomyocytes or progenitor cell types, and instead expressed a group of genes that are specifically expressed in hepatocytes, a cell type stemming from the endoderm layer rather than the mesoderm layer. This small population of cells drove a significant amount of variation in the data (S4 Fig), making it difficult to properly resolve the mesoderm-specific lineages that were the focus of this project. For this reason, the cells assigned to one of the mesoderm-specific lineages (clusters 1–6) were used for downstream analysis (a UMAP embedding of the subset data is shown in S18 Fig).

The log-normalized gene expression was re-centered and re-scaled, and PCA was re-run on specifically these cells to properly focus on the variation among the lineages of interest. The top 3 re-computed PCs were used to calculate a new neighborhood graph, which was used to compute a new embedding to visualize specifically the two cardiac-related differentiating lineages (Fig 2A). The bifurcation into separate cardiac fibroblast and cardiomyocyte lineages can clearly be observed in the PAGA plot (S5 Fig), which was created with the previously described cell type annotations, the re-computed neighborhood graph, and an edge weight threshold of 0.15. This PAGA embedding was used to define the two lineages used for downstream lineage isolation tasks, where all iPSC, mesoderm, cardiac mesoderm, and cardiac progenitor cells are assigned jointly to both lineages, while cardiomyocyte and cardiac fibroblast (terminal cell types) are unique to their corresponding lineage. Finally, four diffusion components were computed from the new neighborhood graph, and diffusion pseudotime was used to assign pseudotime values to cells from both cardiac lineages.

Influence of batch on cell type composition

If a group of cells from the same batch are more likely to display similar cell type composition than cells from different batches, this could indicate that batch effects are confounding cell type annotation. However, since the experimental design places cells from the same cell line and differentiation stage in the same batch, it is important to disentangle similarity in cells due to the same cell line/ differentiation stage (a likely biological effect) from similarity in cells due to batch (a primarily technical effect; see S1 Fig). To do this, we grouped cells by cell line and differentiation day (called a ‘sample’ here), and summarized each group by their cell type proportions. As indicated in S1 Table, each batch contains three different samples representing a different cell line/differentiation day combination. We focused on the impact of batch on terminal cell type composition by analyzing samples from collection day 3 that had been differentiating for 11 or 15 days. At these time points, terminal cell types (cardiomyocyte, CM, and cardiac fibroblast, CF) make up a majority of the cells present. Based on the experimental design, we have 19 pairs of samples such that one sample is at differentiation day 11, the other is at day 15, and the two samples belong to the same batch. If batch effects are causing cells from the same batch to receive similar terminal cell type labels, we would expect these sample pairs to display greater similarity of cell type proportion compared to a background of sample pairs from distinct batches matched for differentiation day. We computed pairwise differences in cell type proportion for both cardiomyocyte and cardiac fibroblast cell types, and compared these to a background (day 11 and day 15 samples from distinct batches) using a t-test. This analysis indicated that samples from the same batch are not more likely to display similar cell type composition than cells from different batches when controlling for cell line and differentiation day (CM p = 0.286, CF p = 0.446).

Pseudobulk expression aggregation and normalization

Although the noisiness of single cell expression profiles necessitates aggregation across cells before dynamic eQTL calling, an improved understanding of the differentiation landscape allows us to pursue an aggregation strategy that mitigates the confounding impact of cellular composition differences and offers greater power than dynamic eQTL calling on bulk samples. Three pseudobulk aggregation schemes were used in this study:

1. Chronological differentiation day binning—This strategy is most directly comparable to bulk RNA-sequencing. Aggregation is performed by taking the sum of SCTransform-corrected counts from all cells from the same differentiation day and individual.
2. Lineage subsetting—Differentiation day binning was performed within each lineage separately. As evidenced by the PAGA graph, all cells up to the progenitor cell type (PROG) are assigned to both lineages, only cells from the terminal cell type (cardiomyocyte or cardiac fibroblast) are unique to one lineage or another.
3. Lineage subsetting & pseudotime binning—After lineage subsetting, cells are partitioned into 16 quantile bins according to pseudotime. We chose 16 bins in order to directly compare to our previous 16 time-point bulk experiment (see S19 Fig). Aggregation then consists of the sum of SCTransform-corrected counts from cells within the same cell line and pseudotime bin.

After pseudobulk aggregation, low-depth samples with library size less than 10,000 were filtered out. In this dataset, this cutoff was equivalent to removing samples with less than or equal to 50 cells. This removed 15 samples from the pseudotime-binned cardiac fibroblast lineage, all corresponding to pseudotime bins 14 or 15 (7 individuals were filtered out of both bins altogether, and one additional individual had only bin 15 filtered out). In the pseudotime-binned cardiomyocyte lineage, this filtered out 5 samples. Remaining samples underwent TMM normalization with singleton pairing through the edgeR package so that expression could be compared across samples for dynamic eQTL calling [69,70]. We then transform the TMM-normalized counts into compute counts per million (CPM) for each sample, and apply log normalization (with the edgeR package, which uses an approach to pseudocount addition that is adapted for library size). These logCPM expression values are used for QTL calling.

Bulk expression normalization

In order to properly compare bulk RNA-seq data to our pseudobulk data, we reprocessed the bulk data from a previous experiment in a way that is intended to most closely match the logcpm pseudobulk expression. For this reason, we used transcripts per million (TPM) instead of previously used reads per kilobase of transcript, per million mapped reads (RPKM). For each sample, we first divided each gene’s counts by the length in kilobase to compute reads per kilobase (RPK), and then fed these adjusted expression values into the same normalization pipeline as was used for pseudobulk counts (which are not biased by gene length)—TMM normalization with singleton pairing and logCPM adjustment, with the edgeR package. Since the input was reads per kilobase rather than counts, this gives logTPM expression values for use in QTL calling.

Sample PCA

To identify primary sources of variation between samples, we ran principal component analysis (PCA) on the gene expression matrix for pseudobulk data. The first principal component is correlated with differentiation time (S20 Fig). For the top 10 PCs, we calculated the percent variance explained of each principal component by each technical factor recorded during sample collection (S14 Fig).

Cell line collapsed PCA

To perform dynamic eQTL calling, we search for changes in gene expression over time that are correlated with a specific genotype. This can be confounded by broad differences between cell lines across the differentiation time course, such as differences in differentiation speed, lineage preference, or technical factors. For example, assume cell lines with genotype G at locus i generally have increasing proportions of cardiomyocytes over time, while cell lines with genotype C at locus i have increased proportions of cardiac fibroblasts over time. In this case, any gene whose expression is upregulated in cardiomyocytes will appear to have a dynamic eQTL at locus i, regardless of any cis-regulatory dynamics related to that gene, which constitute the intended focus of this study.

With single-cell data, we are able to more directly account for some of these factors, namely differentiation speed (with pseudotime binning) and lineage preference (with lineage subsetting). However, it remains useful to control for any broad cell line differences in this more unsupervised fashion, as any broad effects could drive false positive QTL detection.

We used a “cell line collapsed PCA” approach to identify such patterns across the entire time course [23]. To identify cell line collapsed PCs, we rearranged the gene expression matrix from the standard pseudobulk expression quantification such that each row represented expression from one cell line and each column represented a gene at a single time point. After standardizing each column to have zero mean and unit variance, we applied PCA to this matrix to learn a low dimensional representation. Each cell line has a shared loading across all time points, and PCs reflect broad differences in the way cell lines proceed through differentiation. For example, in bulk, it appears that the first cell line PC picks up on differences in differentiation speed between cell lines, while the second cell line PC picks up on differences in terminal cell type preference as defined as the highest total cell type proportion among days 10 to 15 (S9 Fig).

Genotype data

We used previously collected and imputed genotype data for the 19 Yoruba individuals from the HapMap and 1000 Genomes Project [71]. For eQTL analyses, we filtered to variants with no missingness and a minor allele frequency of at least 0.1 across the 19 individuals present.

Dynamic cis-eQTL test selection

We selected which genes to check for dynamic eQTLs based on the following filters:

• Gene must have at least 0.1 CPM in at least 10 bulk/pseudobulk samples
• Gene must have at least 6 counts (reads) in at least 10 samples

Both of these filters were applied separately for each aggregation scheme. We tested all variants within 50kb of the transcription start site of each gene. Transcription start sites were obtained from Gencode’s release 37 [72], basic gene annotation, and matched to mapped genes by Ensembl gene ID. The total number of tests is presented alongside the number of dynamic eQTLs detected in Tables 1 and S2.

Linear dynamic eQTLs using single-cell pseudobulk data

Linear dynamic eQTLs are cis-eQTLs whose effects are linearly modulated by differentiation time. We detected linear dynamic eQTLs with a Gaussian linear model that quantified the interaction between genotype and differentiation time on gene expression, while controlling for the linear effects of both genotype and differentiation time. We also controlled for linear effects of the first five cell line collapsed PCs (see below).

Following the method used in [23], we built a separate linear model for each tested variant-gene pair. Specifically, let t denote the time point (or, for pseudotime binning, the median pseudotime value across cells constituting the pseudobulk sample) of the current sample, c denote the cell line of the current sample, T denote the total number of time points, and C denote the total number of samples. E∈R^CxT denotes the standardized expression matrix for the current gene, G∈R^C denotes the dosage based genotype vector for the current variant, and PC^K∈R^C denotes the Kth cell line collapsed PC vector. We modeled the expression levels as follows:

We used lmFit from the limma package to fit this model, and used a t-test to measure the significance of the genotype and time coefficient (β₁₃).

Bonferroni correction was applied to account for multiple SNPs being tested per gene, and Storey’s q-value was used to control false discovery rates at the gene level, after selecting the most significant dynamic eQTL per gene [73].

Model diagnostics with partial regression and partial residual plots for the example linear dynamic eQTL in Fig 3A are shown in S12 Fig [74].

Nonlinear dynamic eQTLs using single-cell pseudobulk data

To detect dynamic eQTLs whose effect size changes non-linearly with time, we used a second order polynomial basis function over time, which alters the above linear dynamic eQTL model as follows:

Once again, time is either time of collection, or median pseudotime of the sample. As before, we used lmFit from the limma package to fit this model, and this time used a similar t-test to measure the significance of the genotype and quadratic time coefficient (β₂₀). Multiple testing correction was applied as with linear dynamic eQTL calling.

Model diagnostics with partial regression and partial residual plots for the example nonlinear dynamic eQTL in Fig 3C are shown in S12 Fig [74].

Control experiments

We assessed calibration of our dynamic eQTL calling methods with permutations. If we permute the time variable in the interaction term, we do not expect this term to properly capture interactions between genotype and time. For each variant-gene pair, we performed an independent permutation of the time variable in the interaction term, across all (cell line, day) samples. The results of this analysis are shown in S7 Fig.

As another check for confounding factors, we explored the possibility that broad differences between cell lines, such as variation in differentiation speed or trajectory preference, are driving false positive discoveries. If a pair of cell lines share properties such as trajectory preference that confound eQTL analysis, variants where those cell lines share genotype would be more likely to appear as dynamic eQTLs. If this is the case, we would expect pairwise correlation between cell lines according to genotype, across the top 200 dynamic eQTLs, to be higher than expected by chance (compared to a background set of random loci matched for minor allele frequency and distance to transcription start site), as these false dynamic eQTLs are not in fact picking up on distinct cis regulatory patterns but broad cell line patterns. This type of elevated correlation was not observed, suggesting that cell line PCA adequately controls for broad cell line differences (S8 Fig).

We also permuted pseudotime across all cells to investigate sensitivity to uncertainty surrounding pseudotime binning and variation in pseudobulk library sizes. For each of 100 cell-level permutations, we aggregated pseudobulk by the permuted pseudotime quantile bins, and performed dynamic eQTL calling. We used the absolute value of the resulting t-statistics to generate an empirical null distribution for each gene. The empirical p-values generally agreed with the original results up to the significance threshold imposed by the number of permutations performed (p = 0.01, -log10(p) = 2). These results do not suggest inflation of the nominal p-values from the original analysis, and were in fact less conservative than the p-values obtained using Student’s t-distribution (S11 Fig).

Simulations to examine type I errors due to ’double dipping’

We conducted simulations to evaluate potential type I error inflation caused by selective inference. We simulated gene expression data from the following linear mixed model:

Here Y_ijk is the expression of gene k in cell j of individual i, where k = 1,…,1000, j = 1,…,100 and i = 1,…,n. The sample size n is 10 or 20. We assumed one cis-eQTL per gene. To simulate the genotype G_ik, we first generated the minor allele frequency (MAF_k) from Uniform(0.1,0.5) and then generated G_ik~binomial(2, MAF_k).

The other variables included genetic effect size β_k, cell maturity M_ij and its effect size α_k, individual-specific random effect a_ik and error term ϵ_ijk. They were generated from the following distributions:

Note that (a_i1,…,a_i,1000) are i.i.d. across individuals and (ϵ_ij1,…,ϵ_ij,1000) are i.i.d. across individuals and cells, but they are both correlated across genes. To construct a realistic correlation structure, we chose Σ to be the correlation matrix of the expression of 1000 randomly selected genes from our pseudo bulk data. We fixed so that cell maturity and individual specific random effect explained 30% variance of expression and varied . We then generated the genetic effect size β_k~N(0,0.1²) or N(0,0.4²), corresponding to on average 0.4% or 6.3% variance of gene expression explained by genetic effects. The variance of the error term was chosen so that the expression of each gene has unit variance.

We defined the pseudo time in this simulation study to be the first gene expression principal component (PC). We divided the cells into three equal pseudo time bins and averaged expression of the cells for each individual in each pseudo time bin into pseudo bulk expression (). We also calculated the average pseudo time for cells within each pseudo bulk sample, denoted by t_il. We tested two models for dynamic eQTL calling (fitted for each gene k separately): 1) linear mixed model with individual-specific random effects ; 2) linear model without random effects. Type I error was calculated across 1000 genes (S10 Fig). The simulation suggests that a fixed-effect linear model for dynamic eQTL calling, as used in this study, was conservative in the presence of multiple measurements per individual and did not lead to type I error inflation. The more powerful linear mixed model did lead to moderate inflation.

Correlation between bulk and pseudobulk data

We calculated the Pearson correlation of the normalized gene expression matrix from matched bulk RNA-seq data [23] with the normalized gene expression matrix from pseudobulk RNA-seq data. We observed a high correlation of gene expression values between bulk and pseudobulk samples of any given differentiation day (S21 Fig), and a consistent pattern of correlation for all cell lines (S22 Fig).

Bulk dataset deconvolution using single cell data

Cell type deconvolution was performed using CIBERSORTx [53]. The method was first assessed for accuracy using pseudobulk data, where a ground truth is available. Cells from each annotated cell type were split into training (60% of cells) and testing (40%) groups. The annotated Seurat object was subset to training data, and the FindAllMarkers command was used to identify a subset of 404 genes for use in deconvolution. We removed genes that were not measured in bulk, leaving 317 genes for use in deconvolution. A gene expression signature matrix was created from exclusively the training data by taking the sum of SCTransform-corrected counts within each cell type. Normalization of the signature matrix was performed using edgeR: normalization factors were first computed with ‘TMMwsp’ method, then TMM-normalized counts were converted to counts per million. To assess the accuracy of this approach, we then used the same normalization pipeline to aggregate pseudobulk by sample for the testing data, where samples corresponded to a (cell line, differentiation day) combination (S13 Fig). To perform deconvolution of the bulk RNA sequencing data, we used the signature matrix described above and subset the bulk data to the 317 genes contained in the signature matrix.

Cell type interaction eQTLs

To account for variable cell type composition in bulk RNA-seq data, rather than looking for cis-eQTLs whose effects are modulated by time (linear dynamic eQTLs), we looked for those whose effects are modulated by cell type proportion [22]. This mitigates the confounding impact of lineage preference on dynamic eQTL calling, as well as differences in differentiation speed (to the extent that this is captured by cell type proportion). To do so, we replaced the time variable in the dynamic eQTL model with cell type proportion as follows:

Where K_ct is the CIBERSORTx inferred cell type proportions in the sample. Separate models were built for each variant-gene pair, in each cell type except the ‘unknown’ cell type. We additionally explored a model in which we regressed out all cell type proportions (except the unknown cell type, as cell type proportions are constrained to sum to 1).

Note that while all fixed cell type proportion terms are included as covariates, there is only one interaction term for a single cell type proportion. Therefore, once again, separate models were fit for each variant-gene pair, in each cell type except ’unknown’. We found that regressing out additional cell types, not just the one included in the interaction term, led to detection of a greater number of genes with a cell type interaction eQTL (S23 Fig). To check whether these additional covariates were in fact introducing false positive associations between individuals, we measured the pairwise genetic correlation between cell lines among the top hits detected after regressing out additional cell type proportions. We then compared this to the genetic correlation among a set of hits detected before regressing out additional cell type proportions, matched for minor allele frequency. We did not see an increase in genetic correlation among significant tests introduced by incorporation of additional covariates (S24 Fig). However, we did observe a lower replication rate of this expanded set of interaction eQTLs among linear dynamic eQTLs (π₁ = 0.69 and 0.32, respectively, compared to 0.84 and 0.43 under the first model).

We also explored including sample-level principal components as covariates in the linear model:

Where U¹ represents the first sample principal component, as opposed to PC_c¹, the first cell line principal component. Here, we again found that additional covariates led to an increased number of cell type interaction eQTLs detected (S23 Fig): for several cell types (pluripotent cells, mesoderm and progenitor) this figure continued to increase with up to 30 principal components regressed out. With the terminal cell types where more interaction eGenes were detected, the maximum number of hits detected occurred after regressing out 10 principal components. The replication rate among dynamic eQTLs decreased as the number of hits detected increased (π₁ = 0.63 and 0.30 for cardiomyocyte and cardiac fibroblast, respectively, after 5 PCs were regressed out; 0.59 and 0.38 after 10; 0.64 and 0.42 after 20; 0.68 and 0.44 after 30). The results from fitting the first model are reflected in the main text.

Static eQTL calling

The interaction testing frameworks described above are designed to test for significant associations between genotype and some environmental variable (pseudotime, real differentiation time, or cell type proportion). To specifically test for genetic variants associated with gene expression levels across all observed environmental contexts, it is not sufficient to perform a similar hypothesis testing procedure on β₁, the coefficient associated with genotype. This coefficient measures the effect of genotype on expression specifically when the environmental variable is equal to zero, not across all environmental contexts [75]. As a result, estimates of this effect can be highly sensitive to the decision of how to center the environmental variable.

A more reliable way to identify eQTLs that are active across contexts is to perform eQTL calling in each context independently, and then pool information to infer patterns of sharing of regulatory effects across contexts, as done by the mashr package [46]. We performed regular eQTL calling in each cell type, and separately in each day, using the following model:

Where PC_i represents the i^th sample principal component, computed from the single-context (day or cell type) gene expression data. We used both canonical and data-driven covariances, as described in [76], to fit the mash model to a random subset of 250,000 tests. To identify tests that are significant across all contexts, we first subset to the most significant test for each gene in any context. We then used the fitted mash model to count the number of conditions each of these effects is significant in, at a local false sign rate threshold of 0.05. This identified 183 static eQTLs using pseudobulk data aggregated by differentiation day, versus 147 static eQTLs using pseudobulk data aggregated by cell type.

Overlap with published GTEx eQTLs

We used the GTEx v8 release to evaluate replication and overlap of our dynamic eQTLs with variants previously detected in adult tissues. To assess replication in each tissue, we used the qvalue package in R [52] to compute π₁ replication rates among all variant-gene pairs that were declared dynamic eQTLs that were also tested in GTEx. To determine the percentage of variant-gene pairs that were declared both dynamic eQTLs and significant cis eQTLs in GTEx, we incorporated cis eQTLs from all tissues.