Dynamic genetic regulation of gene expression during cellular differentiation

See allHide authors and affiliations

Science  28 Jun 2019:
Vol. 364, Issue 6447, pp. 1287-1290
DOI: 10.1126/science.aaw0040

Variants affect gene expression over time

Genetic variation drives the spectrum of human phenotypes but in some cases has been linked to pathological conditions. Strober et al. set out to explore how genetic diversity regulates gene expression in cell differentiation over time. They examined expression quantitative trait loci (eQTLS)—genetic variants correlated with gene expression—in induced pluripotent stem cell lines from 19 Yoruban individuals at 16 time points during differentiation into cardiomyocytes. They identified hundreds of developmental stage–specific eQTLs. Thus, the impact of genetic variations depends on the developmental context or tissue under observation.

Science, this issue p. 1287


Genetic regulation of gene expression is dynamic, as transcription can change during cell differentiation and across cell types. We mapped expression quantitative trait loci (eQTLs) throughout differentiation to elucidate the dynamics of genetic effects on cell type–specific gene expression. We generated time-series RNA sequencing data, capturing 16 time points during the differentiation of induced pluripotent stem cells to cardiomyocytes, in 19 human cell lines. We identified hundreds of dynamic eQTLs that change over time, with enrichment in enhancers of relevant cell types. We also found nonlinear dynamic eQTLs, which affect only intermediate stages of differentiation and cannot be found by using data from mature tissues. These fleeting genetic associations with gene regulation may explain some of the components of complex traits and disease. We highlight one example of a nonlinear eQTL that is associated with body mass index.

Genetic variants that alter gene regulation play an essential role in the genetics of human disease and other complex phenotypes (1, 2). Large studies have identified thousands of genetic loci associated with complex diseases, most of which are in noncoding regions of the genome and therefore are putatively involved in gene regulation (2). Expression quantitative trait locus (eQTL) analysis has shown that many disease-associated loci influence the regulation of nearby genes (3, 4), but a substantial fraction of disease-associated loci still remain unexplained (5, 6).

Much effort has been dedicated to mapping and identifying eQTLs across tissues and cell types, as the regulatory impact of disease-associated loci may be most evident in cell types relevant to each disease. Regulatory genetic effects can also be time point specific or environment dependent (7, 8) and may influence temporal programs of gene regulation. Yet almost all studies of the genetics of gene regulation, including the multitissue Genotype-Tissue Expression (GTEx) project (7), involve data collected at a single time point, usually from adult individuals. Dynamic gene expression data can add another dimension to eQTL analysis, allowing identification of genetic variants with transient effects that may not have been found in analysis of static data.

We took advantage of a panel of induced pluripotent stem cell (iPSC) lines from 19 individuals to investigate high-resolution temporal genetic effects on gene regulation over time during cardiomyocyte differentiation. Specifically, we collected gene expression data throughout the differentiation from iPSCs to cardiomyocytes in 19 well-characterized human Yoruba HapMap cell lines (9). For each cell line, RNA was extracted and sequenced every 24 hours for 16 days to capture the entire differentiation process; in total, we sequenced 297 RNA samples (figs. S1 and S2). Combined with available whole-genome sequences and genotype data for each cell line, these data provide a resource with which to investigate how gene expression and genetic regulation change throughout cardiomyocyte differentiation with high temporal resolution.

During iPSC culturing, differentiation, RNA extraction, and processing for sequencing, we recorded extensive metadata on each sample (table S1). Quality controls and filtering yielded 16,319 genes for downstream analysis (10). After standardization and normalization of the RNA sequencing (RNA-seq) data (10), we evaluated the contribution of potential confounders to overall variation in our data, confirming that our study design was effective (fig. S3). We also used replicates from an independent differentiation to confirm that the gene expression patterns we observed in our iPSCs and iPSC-derived cardiomyocytes are robust with respect to variance that may be associated with the differentiation procedure (fig. S4) (9, 10).

We evaluated the efficiency of our differentiation by fluorescence-activated cell sorting (table S2) and by considering the time-course expression of known cell type–specific marker genes (11, 12) (fig. S5). As expected (12), cardiomyocyte purity and the expression of lineage marker genes are variable across our samples. This variability between cell lines was observed across the entire time course, although the effect of differentiation time is the primary source of variation in the data (Fig. 1A and figs. S3 and S6).

Fig. 1 Gene expression trends throughout cardiomyocyte differentiation.

(A) The first two gene expression principal component (PC) loadings for all 297 RNA-seq samples across cell lines, where each sample is colored according to the day of collection. (B) Predicted cell line cluster expression trajectories for 20 gene clusters according to split-GPM. Many gene clusters (8, 11, 15, 16, and 20) exhibit periodic expression trajectories that correspond with cell culture media changes.

We characterized global patterns of gene expression across time by applying split-GPM, an unsupervised probabilistic model that infers time-course trajectories of gene expression using Gaussian processes, while simultaneously performing clustering of genes and cell lines (10). Using this approach, we identified two clusters of cell lines that displayed broad differences in the expression patterns of multiple clusters of genes; in each gene cluster, genes exhibit shared expression changes over time. The assignment of cell lines to clusters is robust with respect to the parameters we tested, such as the number of inferred gene clusters (fig. S7).

The two cell line clusters we identified differ in the efficiency of cardiomyocyte differentiation. Cell lines in the first (larger) cluster display greater troponin expression levels in the final six time points of differentiation (P = 0.014, Wilcoxon rank-sum test). The expression of a group of genes enriched for myogenesis also increases by a greater magnitude over time in cell lines in the first cluster (Bonferroni P = 9.29 × 10−14) (gene cluster 2 in Fig. 1B) (13). Cell lines in the second, smaller cluster show high expression of genes related to KRAS activation (Bonferroni P = 0.005; gene cluster 4 in Fig. 1B), which is associated with increased self-renewal of undifferentiated iPSCs and decreased neuronal differentiation propensity (14). Other gene clusters illuminate broad changes in gene expression over time, such as a transient rise in MYC and E2F target genes in the early days of differentiation (gene cluster 13 in Fig. 1B; table S3). Together, this analysis documents patterns of gene expression trajectories over time and captures differences among our cell lines that are not obvious from the individual time point data alone.

Next, we evaluated the impact of genetic variation on gene regulation in our system. We used WASP software (15) to identify cis-eQTLs in the data from each time point independently (10). To control for latent confounders in the independent analysis of data from each time point, we included the first three expression PCs using data from samples of the corresponding time point as covariates (figs. S8 and S9, A and B). At an empirical false discovery rate (eFDR) of 5%, we identified a median of 111 genes (range: 71 to 231) with at least one eQTL in each time point (figs. S9C and S10). As expected, the eQTLs we identified early in the time course replicated in data from iPSCs, whereas eQTLs from later time points were better supported by data from iPSC-derived cardiomyocytes (both P < 0.001, linear regression) (Fig. 2A) (9).

Fig. 2 eQTL patterns during cardiomyocyte differentiation.

We limit this analysis to genes with at least one significant eQTL (WASP combined haplotype test; eFDR ≤ 0.05) across time points. If a gene has more than one significant eQTL, we select a single variant for that gene with the smallest geometric mean P value across all 16 time points. (A) Spearman correlation of P values between eQTLs from each day (x axis) and existing iPSC (gray) and iPSC-derived cardiomyocyte (red) eQTLs. (B) Spearman correlation of eQTL P values for each pair of days. (C) Factors identified via sparse matrix factorization of eQTL-log10 P values using three latent factors and an L1 penalty of 0.5.

We computed the correlation of the significant eQTL summary statistics for each pair of time points (Fig. 2B). We observed that correlation between eQTL summary statistics increases as the distance between time points decreases (P ≤ 2 × 10−16, linear regression). Although this observation is intuitive, it indicates that the dynamic impact of genetic variation on gene regulation in our data is not random and is related to the temporal process of cardiomyocyte differentiation.

To more formally quantify the temporal structure of genetic regulation throughout differentiation, we performed sparse non-negative matrix factorization on the matrix of significant eQTL summary statistics from all time points (10). The learned factors capture genetic signal that is largely specific to a subset of differentiation time (Fig. 2C), a pattern that is robust with respect to the number of latent factors or sparse prior choice (fig. S11).

Our analysis indicates that temporal structure dominates the patterns of genetic association with gene expression in our data. However, the observation that most significant nondynamic eQTLs can be identified in only a few time points (median of 2) (fig. S12) is most likely explained by incomplete power to identify eQTLs in each time point independently. To robustly identify dynamic eQTLs whose effect varies significantly over time, leveraging power across all time points (Fig. 3A), we used a Gaussian linear model applied jointly to data from the entire experiment. Specifically, we quantified the effect of interactions between genotype and differentiation time on gene expression, controlling for linear effects of both differentiation time and genotype. In addition, we accounted for the systematic differences in differentiation trajectories identified between cell lines (Fig. 1B, figs. S13 to S16, and table S4) (10), which would otherwise lead to false positives in our analysis. Using this approach, we identified 550 genes with a significant dynamic eQTL (eFDR ≤ 0.05) (figs. S17 to S20 and table S5).

Fig. 3 Dynamic eQTLs detect genetic regulatory changes caused by cardiomyocyte differentiation.

(A) Linear interaction association between genotype (color) of rs11124033 and time point (x axis) on residual gene expression (cell line effects regressed on expression) of FHL2 (y axis). (B) Enrichment of dynamic eQTLs in cell type–specific chromHMM enhancer elements relative to 1000 sets of randomly selected matched-background variants. Dynamic eQTLs were classified as either early or late. (C) Nonlinear interaction association between genotype (color) of rs28818910 and time point (x axis) on residual gene expression of C15orf39 (y axis). (D) Nonlinear interaction association significance of all variants tested within 50 kb of the C15orf39 transcription start site with expression of C15orf39 (green) and GWAS significance for BMI of variants in the same window (blue). Vertical line depicts genomic location of the most significant nonlinear dynamic eQTL (rs28818910) for C15orf39.

We classified the 550 dynamic eQTL as early (eQTL effect size decreasing over time), late (eQTL effect size increasing over time), or switch (eQTL effect size exhibiting different directions of effect over time) (fig. S21) (10). We found that the early dynamic eQTLs are enriched for chromHMM enhancer elements annotated in iPSC Roadmap Epigenomics cell types but not in heart-related cell types (16, 17). In turn, late dynamic eQTLs are enriched for chromHMM enhancer elements annotated in heart-related Roadmap Epigenomics cell types but not in iPSCs (Fig. 3B and fig. S22). These observations indicate that dynamic eQTL mapping can capture temporal changes in cellular gene regulation reflecting changes in regulatory element activity as the cell cultures differentiate.

The observation that we are able to capture the function of cell type–specific regulatory elements prompted us to consider dynamic eQTLs in other contexts. We found that dynamic eQTLs are enriched for genes with roles in myogenesis (Bonferroni P = 0.0019, Fisher’s exact) (table S6) (13) and also show significant enrichment for genes related to dilated cardiomyopathy (P = 0.001, Fisher’s exact) (table S7) (10, 18). Two significant dynamic eQTLs in particular, rs7633988 and rs6599234 (in strong linkage disequilibrium; coefficient of determination, R2 = 0.93), are genome-wide association study variants for QRS duration and QT interval, respectively (fig. S23) (19, 20). Both variants show an association with the expression levels of SCN5A, which is involved in the creation of sodium channels and is in the dilated cardiomyopathy gene set (21). Another dynamic eQTL, rs11124033, associated with the expression of FHL2 (Fig. 3A), is also associated with dilated cardiomyopathy. This variant lies in a Roadmap Epigenomics chromHMM promoter element annotated in heart-related cell types but not in iPSCs (16, 17). None of these examples were identified as eQTLs in the nondynamic QTL analysis of each time point from our dataset or in the GTEx heart tissue data (7).

Finally, we sought to identify a wider range of dynamic regulatory patterns, including nonlinear associations, such as when a genetic effect increases in magnitude in the middle of the time course before decreasing or disappearing. To identify nonlinear dynamic eQTLs, we expanded our linear model using a second-order polynomial basis function (10). We acknowledge that our study is underpowered to expand to a more general class of nonlinear dynamic eQTLs that do not assume a continuous effect of differentiation time (fig. S24) (10).

We identified 693 genes with a nonlinear dynamic eQTL (eFDR ≤ 0.05) (figs. S17B and S19B and table S8), 28 of which have their strongest genetic effect in the middle of the differentiation time course (middle dynamic eQTLs) (fig. S25) (10). Twenty-five of these middle dynamic eQTL genes and their strongest associated variant are not identified as eQTLs in our nondynamic QTL analysis in either iPSCs (day 0) or cardiomyocytes (day 15).

In one example of a nonlinear dynamic eQTL, rs8107849 is associated with the expression of ZNF606 with a larger magnitude of effect during days 4 through 11 (fig. S26). The rs8107849 locus does not lie in iPSC or heart-related chromHMM regulatory regions and was not identified in our analysis as a nondynamic eQTL at any time point. Although ZNF606 is known to have a role in differentiation of chondrocytes (22), it is possible this is a conserved process involved in the differentiation of additional cell types, including cardiomyocytes. Another nonlinear dynamic eQTL reveals an association between rs28818910 and C15orf39. The rs28818910 variant is also associated with body mass index (BMI) (P < 6.07 × 10−9, reported) (Fig. 3, C and D) (23) and weakly associated with red blood cell count (P < 1.48 × 10−6, reported) (24). This dynamic eQTL and both traits show similar patterns of association across the region (fig. S27). The rs28818910 locus is associated with interindividual differences in gene expression only during intermediate stages of differentiation; it does not lie in annotated regulatory elements of either iPSCs or cardiomyocytes and is not identified as an eQTL in iPSCs, mature cardiomyocytes, or either of the two GTEx heart tissues. Thus, this is an example of a temporary dynamic regulatory effect that may have phenotypic consequences.

Our time-course study design allowed us to identify hundreds of dynamic eQTLs throughout the differentiation of human iPSCs to cardiomyocytes. Dynamic eQTLs, in particular those with nonlinear effects, may often be transient and will not be found in studies that only consider gene expression data from either stem cells or mature tissues and cell types. Many of our dynamic eQTLs lie in regions without known regulatory annotations, as functional studies have focused on static cell types. Thus, these loci may have previously unknown regulatory effects, which could be followed up with further functional validation in relevant intermediate time points. The dynamic genetic effects identified in our study, or in future time-series genomic datasets, will provide a resource for investigating mechanisms underlying disease associations that cannot be characterized based on studies of terminal cell types.

Supplementary Materialss

Materials and Methods

Figs. S1 to S27

Tables S1 to S8

References (2834)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: We thank H. Kyung Im for help with GWAS analysis and J. K. Pritchard for providing comments on the manuscript. Funding: Y.G. and A.B. were supported by NIH/NIGMS R01GM120167. R.E. was supported by the NIH MSTP Training Grant T32GM007281. K.R. was supported by NIH GRTG 5T32GM007197 and AHA Predoctoral Fellowship 18PRE34030197. The computational resources were provided by the University of Chicago Research Computing Center. Author contributions: Y.G. and A.B. conceived the study. R.E. and K.R. performed the experiments. K.T. developed split-GPM. B.J.S. analyzed the data, with assistance from N.K. All authors wrote the paper. Y.G. and A.B. supervised this project. Competing interests: The authors declare no competing interests. Data and materials availability: The fastq files and RNA-seq counts are available at GEO under accession GSE122380. Files containing nondynamic and dynamic eQTL summary statistics are available at (25). The split-GLM package can be found at (26). Scripts used for this analysis can be found at (27).

Stay Connected to Science

Navigate This Article