Circadian networks in human embryonic stem cell-derived cardiomyocytes.

Cell-autonomous circadian oscillations strongly influence tissue physiology and pathophysiology of peripheral organs including the heart, in which the circadian clock is known to determine cardiac metabolism and the outcome of for instance ischemic stress. Human pluripotent stem cells represent a powerful tool to study developmental processes in vitro, but the extent to which human embryonic stem (ES) cell-derived cardiomyocytes establish circadian rhythmicity in the absence of a systemic context is unknown. Here we demonstrate that while undifferentiated human ES cells do not possess an intrinsic functional clock, oscillatory expression of known core clock genes emerges spontaneously during directed cardiac differentiation. We identify a set of clock-controlled output genes that contain an oscillatory network of stress-related transcripts. Furthermore, we demonstrate that this network results in a time-dependent functional response to doxorubicin, a frequently used anti-cancer drug with known cardiotoxic side effects. Taken together, our data provide a framework from which the effect of oscillatory gene expression on cardiomyocyte physiology can be modeled in vitro, and demonstrate the influence of a functional clock on experimental outcome.

Introduction

The circadian clock is a conserved time‐keeping system that regulates numerous body features such as behavior, metabolism, body temperature, tissue regeneration, and organ homeostasis in a diurnal manner 1. In the heart, the role of 24‐h rhythmicity is illustrated by oscillations in heart rate, blood pressure, and cardiac output 2, 3, 4, 5, 6. The circadian clock comprises a central clock in the brain (the suprachiasmatic nucleus (SCN)) and peripheral clocks that are present in almost all organs. The SCN is mainly entrained by light and subsequently synchronizes the peripheral clocks via neural and humoral factors 7. Interestingly, peripheral clocks function in a cell‐autonomous manner. When ablating the SCN, these clocks remain functional and even synchronized when subjected to a restricted feeding regime 8, 9. Autonomous rhythmicity is underscored by the persistence of circadian rhythms in in vitro cultured cells.

The molecular mechanism that underlies the core clock machinery consists of a transcriptional/translational feedback loop in which a heterodimer of BMAL1 and CLOCK drives rhythmic transcription of downstream genes. These include other core clock genes (period 1 (PER1), PER2, PER3, cryptochrome 1 (CRY1), CRY2, RORα/β, REV‐ERBα/β) as well as clock‐controlled genes (CCGs) that determine circadian organ physiology in a tissue‐specific manner. In the murine heart, ~6–12% of the expressed genes have a circadian expression pattern 10, 11, 12, 13, 14. Rhythmicity is essential for human tissue homeostasis as well, as highlighted by the fact that genetic or environmental (e.g., shift‐work) perturbation of the circadian clock results in a vast array of malignancies such as sleep disorders, inflammation, cancer 15, impairment of regenerative capacity 16, 17, metabolic disorders 18, 19, 20, and cardiovascular diseases 13, 21, 22, 23, 24. In addition, the onset of multiple malicious cardiac events is known to follow a diurnal pattern. Myocardial infarction 25, 26, arrhythmias 27, 28, and sudden cardiac death 29, 30 show a higher incidence in the sleep‐to‐wake transition in humans. The important role of circadian rhythmicity in cardiac injury and regeneration is further solidified by genetic experiments in mice in which a cardiomyocyte‐specific mutation of the Clock gene has been shown to blunt the heart's response to induced ischemic damage 31. Accordingly, clinical studies revealed that infarcts were larger and led to increased reduction in cardiac function when occurring in the sleep‐to‐wake transition 32, 33, 34.

Human pluripotent stem cell‐derived cardiomyocytes have emerged as a potential cellular source for replacement therapies. In addition, human ES cell‐derived as well as induced pluripotent stem cell‐derived cardiomyocytes are increasingly used for disease modeling and drug testing 35. While circadian rhythms play an essential role in cardiomyocyte function in vivo, nothing is known about circadian control of gene expression in pluripotent stem cell‐derived cardiomyocytes, which are often used to model cardiac function and disease.

Here we analyze temporal clock gene expression networks in human ES cells and ES cell‐derived cardiomyocytes. We demonstrate that circadian rhythmicity is absent in human ES cells and is established progressively during directed cardiac differentiation. The identified oscillatory networks are shown to significantly influence the function of human ES cell‐derived cardiomyocytes and determine their response to externally applied stressors. Our findings underscore that circadian rhythmicity can affect experimental outcome, which may have important ramifications for processes such as timed cell‐based therapy.

Results

Human embryonic stem cells express clock genes in a non‐oscillatory manner

Nearly all cells in the human body possess a functional clock as indicated by circadian rhythmicity of core clock gene expression. However, whether human embryonic stem (ES) cells display a functional circadian clock is unknown. Therefore, we compared global expression levels of six core clock genes ARNTL (coding for and henceforth referred to as BMAL1), PER2, CRY1, CRY2, CLOCK, and NR1D1 between pluripotent human ES cells and differentiated human osteosarcoma U2OS cells, a cell line known to possess a functional clock as confirmed here by anti‐phasic Bmal1‐dLuc and Per2‐dLuc signal after transduction of promoter‐based lentiviral luciferase constructs 36, 37 (Fig 1A and B). Transcripts of all genes could be detected in both cell types, with five out of six core clock genes having higher expression levels in U2OS cells compared to human ES cells (Student's t‐test, P < 0.05; Fig 1C). Only CRY1 showed a trend toward higher expression in human ES cells (Student's t‐test, P = 0.0506; Fig 1C), which is in line with previously reported observations in mouse ES and NIH3T3 cells 38. Protein levels of BMAL1, CRY1, and CLOCK were detected by Western blot (Fig 1D, left and Fig EV1) at corresponding levels to their mRNA transcripts (Fig 1D, right). From this, we conclude that while core clock genes are expressed and translated into proteins in human ES cells, this occurs with a different stoichiometry in comparison with U2OS cells and most likely also other differentiated cell types.

Figure 1

Non‐oscillatory expression of clock genes in human ES cells

Raw lentiviral promoter‐based luciferase reporter bioluminescence in U2OS cells after dexamethasone synchronization. Bioluminescence was measured with a LumiCycle32. Values are relative to T0.

Detrended bioluminescent signals measured in (A).

BMAL1, PER2, CRY1, CRY2, CLOCK, and NR1D1 expression levels in human ES cells compared to U2OS cells as determined by qRT–PCR. Expression levels were normalized to PPIA and compared between cell types using an unpaired two‐tailed Student's t‐test (ns: not significant, *P < 0.05, **P < 0.005, and ***P < 0.0005). Data are represented as mean ± s.e.m. of three independent replicates.

Western blot for BMAL1, CRY1, and CLOCK. Protein levels were quantified and normalized to β‐ACTIN.

qRT–PCR analysis of BMAL1 and PER2 expression over 48 h at a 4‐h interval in human ES cells. Circadian oscillations were analyzed using the RAIN algorithm, and the significance of rhythmicity across 48 h is indicated (ns: not significant). Data are represented as mean ± s.e.m. of three independent replicates.

Bmal1‐dLuc and Per2‐dLuc values in synchronized human ES cells across 76 h measured by LumiCycle32. Representative tracks are shown. Values are relative to T0.

Figure EV1

Protein levels of clock genes in human ES cells versus U2OS cells

Original Western blots for BMAL1, CRY1, and CLOCK protein.

The presence of core clock proteins in human ES cells led us to investigate their possible rhythmic expression over time. To this end, human ES cells were synchronized with forskolin 39, and BMAL1 and PER2 mRNA levels were measured every 4 h over a period of 48 h using qRT–PCR. Significance of 24‐h rhythmicity was assessed using RAIN, a nonparametric method detecting arbitrary wave forms in biological data 40. No apparent oscillatory expression pattern could be identified over the course of 2 days (RAIN, BMAL1, P > 0.99 and PER2, P = 0.97; Fig 1E). Additionally, to assess BMAL1 and PER2 transcription over time in human ES cells, we transduced human ES cells with Bmal1‐ and Per2‐promoter‐based lentiviral luciferase constructs 36, 37. After synchronization, no rhythmic bioluminescence was observed (Fig 1F). Therefore, clock genes are expressed in human ES cells, but in a non‐circadian manner.

Human embryonic stem cell differentiation toward cardiomyocytes

Multi‐lineage differentiation of human ES cells has proven extremely valuable to understand developmental processes as well as to provide clinically relevant populations for cell‐based therapy and drug testing 41. To assess the presence of a functional clock upon differentiation, circadian rhythmicity was analyzed at two additional stages (D15 and D30) during directed differentiation of human ES cells toward cardiomyocytes (Fig 2A). To allow for the identification of early cardiac cells, we made use of a NKX2‐5‐eGFP (Homeobox protein NKX2‐5‐eGFP) reporter human ES cell line 42. Cardiac differentiation of human ES cells typically yields significant contribution of cardiomyocytes to the total population of cells 43, 44, 45, which was also seen here with ~50% cardiomyocytes around D15 as defined by FACS for cTNNT2 (cardiac Troponin T) (Fig EV2A). Different stages were characterized by clear changes in marker gene expression (Fig 2A). At day 0 (D0), cells expressed the pluripotency markers NANOG and POU5F1, both at the mRNA and protein level (Figs 2A and B, and EV2B). Upon differentiation, pluripotency factors quickly decreased and the expression of cardiac markers, such as NKX2‐5 and ACTC1 (actin alpha cardiac muscle 1), was observed in both (D15 and D30) spontaneously beating cultures as measured by qRT–PCR (Fig 2A). Immunostainings for α‐sarcomeric ACTIN and cTNNT2 confirmed sarcomeric structures at D15 and D30 (Fig 2B). Additionally, staining for MEF2C (Myocyte Enhancer Factor 2C) and GFP, to assess the presence of NKX2‐5‐eGFP‐positive cells, revealed the abundance of cardiomyocytes at D15 and D30. While the early cardiomyocyte progenitor marker ISL1 (insulin gene enhancer protein ISL‐1) was highly expressed at D15, maturation markers such as KCNJ2 (inward rectifier potassium channel 2) and SERCA2A (sarcoplasmic/endoplasmic reticulum Ca2+ ATPase) 44, 46 were the highest at D30 (Fig 2A). These results validate the in vitro transcriptomic maturation of these cells between D15 and D30, and confirm that the different stages represent distinct cardiac states that can be used to assess the presence of a functional clock across the transition from human ES cells to cardiomyocytes.

Figure 2

Characterization of distinct stages of cardiac differentiation

Schematic of the directed cardiac differentiation and the three different stages used in this study. hESCs: human embryonic stem cells. CM: cardiomyocyte. Below are mRNA expression levels of pluripotency and cardiac markers as measured by qRT–PCR. Expression levels were normalized to a non‐oscillatory housekeeping gene (PPIA). Data are represented as mean ± s.e.m. of three independent replicates.

Immunostaining for pluripotency markers OCT4 and NANOG (yellow) in human ES cells. cTNNT2 and sarcomeric actin (red) stainings reveal clear sarcomeric structures at all cardiac stages. Cardiomyocyte nuclei were stained for MEF2C (red), NKX2‐5‐eGFP‐positive cells were stained with anti‐GFP (green) and nuclei with Hoechst (blue).

Figure EV2

Cardiomyocyte quantification and characterization by FACS and qRT–PCR

FACS staining for cTNNT2 in D10 and D15 human ES cell‐derived cardiac cells. IgG1 isotype was used as a control.

POU5F1 mRNA expression levels in human ES cells and D15 as well as D30 cardiac cells, as measured by qRT–PCR. Expression levels were normalized to a non‐oscillatory housekeeping gene (PPIA). Data are represented as mean ± s.e.m. of three independent replicates.

Rhythmic expression of clock genes emerges during cardiac differentiation

Whether and when human ES cells develop a functional clock upon differentiation, in the absence of systemic cues, is unknown. To investigate the possibility of an emerging clock, we compared mRNA levels of BMAL1, PER2, and CLOCK in D0 human ES cells and D15 and D30 human ES cell‐derived cardiac cells. Even though BMAL1, PER2, and CLOCK were expressed at all stages, their expression increased significantly from D0 to D15 and/or D30 (ANOVA with Bonferroni correction, P < 0.05; Fig 3A). This indicates that core clock gene expression levels gradually increase during directed cardiac differentiation.

Figure 3

Rhythmic clock gene expression emerges in (matured) cardiac cultures

Data information: Measurements in (F–H) were performed with a LV200 microscope.

BMAL1, PER2, and CLOCK expression levels at D0, D15, and D30 during directed cardiac differentiation as determined using qRT–PCR. Data are represented as mean ± s.e.m. of three independent replicates. Significant expression differences were tested by one‐way ANOVA, followed by a Bonferroni post hoc test (ns: not significant, *P < 0.05, **P < 0.005, ***P < 0.0005).

qRT–PCR analysis of BMAL1 and PER2 expression over 48 h at a 4‐h interval in cardiac cells at D15. Expression levels in (A) and (B) were normalized to PPIA. Data are represented as mean ± s.e.m. of three independent replicates. Significance of rhythmicity across 48 h was analyzed using the RAIN algorithm and is indicated (ns: not significant, ***P < 0.0005).

Promoter‐based destabilized luciferase (dLuc) reporter assay of the Bmal1 and Per2 promoter in synchronized cardiac cells at D15. Values are relative to T0. Measurements were performed using a LumiCycle32.

Similar analysis as in (C) for D30.

Detrended Bmal1‐dLuc and Per2‐dLuc luciferase signal measured in (D).

Detrended Bmal1‐dLuc and Per2‐dLuc bioluminescence in NKX2‐5‐eGFP+ sorted and synchronized human ES cell‐derived cardiomyocytes at D30.

Single‐cell analysis of Per2‐dLuc bioluminescence in sorted eGFP‐positive and synchronized D30 human ES cell‐derived cardiomyocytes.

Representative Per2‐dLuc signal in five single D30 human ES cell‐derived cardiomyocytes over the course of 48 h.

To assess rhythmicity of clock gene expression at D15 and D30, cells were synchronized using dexamethasone 47 and three independent RNA samples were collected every 4 h over a period of 48 h (Fig EV3). BMAL1 and PER2 levels were analyzed by qRT–PCR to determine whether their expression oscillated in an anti‐phasic manner, a hallmark of a functional molecular circadian clock. Similar to undifferentiated human ES cells (Fig 1E), no clear circadian pattern was observed in the early stage D15 human ES cell‐derived cardiomyocytes (D15; RAIN, P = 0.095, P = 0.68 for BMAL1 and PER2, respectively; Fig 3B). Matured cardiac cells, however, showed significant oscillations for PER2 but not BMAL1 (D30; RAIN, P = 1.4E‐8 and P = 0.81; Fig 3B). To further validate the emergence of a functional clock, cultures were transduced with Bmal1‐dLuc and Per2‐dLuc lentiviral reporters. At D15, after synchronization, a small induction of oscillatory Per2‐based luciferase signal could be detected (Fig 3C), which is in line with previously described observations of Per2 as an early oscillator upon retinoic acid induced differentiation in mouse ES cells 48. In D30 synchronized populations, typical anti‐phasic oscillatory Per2‐ and Bmal1‐driven bioluminescence levels were observed, which confirms the presence of a clock at D30 (Figs 3D and E, and EV4A and B).

Figure EV3

Setup of CEL‐Seq experiment

Schematic of RNA samples that were processed for qRT–PCR and/or CEL‐Seq. ZT: Zeitgeber and RPM: reads per million.

Figure EV4

Human ES cell‐derived cardiomyocytes show circadian oscillation in Bmal1‐ and Per2‐dLuc

Promoter‐based destabilized luciferase (dLuc) reporter assay of the Bmal1 (red) and Per2 (black) promoter in synchronized human ES cell‐derived cardiac cells at D30 (replicate for Fig 3D). Values are relative to T0. Right: Detrended Bmal1‐dLuc and Per2‐dLuc luciferase signal. Measurements were performed using a LumiCycle32.

Similar analysis as in (A) on an additional independent replicate.

Brightfield and fluorescent image of D30 GFP+ sorted human ES cell‐derived cardiomyocytes used for analysis in Fig 3F (left/red).

FACS plot, brightfield, and fluorescent image of D30 GFP+ sorted human ES cell‐derived cardiomyocytes used for single‐cell bioluminescent analysis of Per2‐dLuc signal in Fig 3F (right/black) and Fig 3G and H.

Similar analysis as in (A) for D45 human ES cell‐derived cardiac cells.

In order to verify the contribution of cardiomyocytes to the observed oscillatory pattern in our cardiac cultures, NKX2‐5‐eGFP+ cells were purified via FACS. After sorting, strong Bmal1‐dLuc and Per2‐dLuc rhythmicity was detected in synchronized human ES cell‐derived cardiomyocytes (Figs 3F and EV4C and D). In addition to these observations in a sorted population, circadian Per2‐dLuc patterns were also found in bioluminescence recordings of single NKX2‐5‐eGFP+ human ES‐derived cardiomyocytes (Figs 3G and H, and EV4D), which further confirms the presence of a functional clock in D30 cardiomyocytes. To question whether circadian rhythmicity would persist during culture, 45 days old cardiac cultures were analyzed using real‐time reporter‐based luciferase measurements. Significant Per2‐ and Bmal1‐dLuc oscillations were observed (Fig EV4E). These results indicate that human ES cells develop a functional clock upon directed cardiac differentiation, with robust oscillations at D30 that persist in older in vitro cultures.

Human ES cell‐derived cardiomyocytes show a network of stress‐related clock output genes

A functional circadian clock translates into the oscillatory expression of clock‐controlled genes (CCGs). Gene expression profiling in numerous cell types and tissues has shown that around 3–16% of the transcriptome exhibits circadian rhythmicity 14. Conform differing physiological demands of organs, oscillating output genes vary per tissue. For the murine heart, ~6–12% of the expressed genes were shown to oscillate in a 24‐h manner 10, 11, 12, 13, 14. To identify CCGs during in vitro cardiomyocyte differentiation, genome‐wide mRNA levels were assessed by mRNA sequencing of purified RNA using CEL‐Seq, a previously described RNA profiling technique based on sequencing the 3′UTR of mRNAs, generating one read per transcript 49. We first compared the overall transcriptional profile of matured cardiac cells (D30) to that of human ES cells and D15 cultures 48 h post‐synchronization (Fig EV5A and Table EV1). Lowly expressed genes are typically not picked up robustly when analyzing highly multiplexed CEL‐Seq data at relatively low sequencing depth. To control for this, genes with an average of < 3 RPM (reads per million) across all time points were not used for further analysis. Based on ~14,000 genes with an expression level of more than 3 RPM at one of the stages, Spearman's rank correlation coefficients (ρ) showed that transcriptional programs were substantially different between cardiac cells at D15 and D30 (ρ = 0.53; Fig EV5A–C). Observed changes between states were consistent between our qRT–PCR and CEL‐Seq analyses, as indicated for several marker genes, which highlights the reliability of our sequencing datasets (Figs 2A and EV5D). Indeed, increased MYH7/MYH6 levels and multiple other markers (e.g., MYL2, PLN, and KCNJ2) confirm in vitro cardiomyocyte maturation as well as generally higher clock gene expression across differentiation (Fig EV5C and D).

Figure EV5

CEL‐Seq‐based characterization of transcriptome maturation and association between read coverage and oscillatory expression

Spearman's rank correlation with average linkage based on average RPM‐normalized CEL‐Seq read counts of three ZT48 replicates. Genes (n = 13,770) with an average ZT48 expression level of more than 3 RPM (reads per million) at one stage (either D0, D15, or D30) were used for analysis. Spearman's rank correlation coefficients are indicated in white.

Heatmap of average ZT48 CEL‐Seq read counts for 13,770 genes used in (A). Genes were sorted according to human ES cell expression from high (top) to low (bottom). Log2 RPM values were plotted.

MYH7/MYH6 CEL‐Seq read count ratio in D15 and D30 human ES cell‐derived cardiac cultures (average of three ZT48 replicates) as a measure for maturation.

Heatmaps of average ZT48 CEL‐Seq read counts in human ES cells (D0) and human ES cell‐derived cardiac cells (D15 and D30). The left heatmap shows pluripotency markers, the two middle heatmaps depict cardiomyocyte and cardiomyocyte maturation markers, and the right panel comprises core circadian clock genes. Log2 RPM values were plotted.

Boxplots depict average log2 RPM‐normalized CEL‐Seq counts across all time points for D15 and D30 oscillatory and non‐oscillatory genes. 10,089 genes with average RPM values above 3 across all time points (ZT4‐ZT48) on which JTK analysis was done in both D15 and D30 were used. Bottom and top of the boxes are the first and third quartiles. The line within the boxes represents the median and whiskers denote the interval within 1.5? the interquartile range from the median. Outliers are plotted as dots.

Boxplots depict average log2 RPM‐normalized CEL‐Seq counts across all time points (ZT4‐ZT48) for D15‐specific, D30‐specific, and shared cardiac oscillators. Boxplot representation is as described in (E).

To assess the possible presence of oscillatory transcripts at D15 and the identity of CCGs in D30 cardiac cells, in which a functional clock was found (Fig 3), three independent RNA samples were collected every 4 h over a period of 48 h and sequenced using CEL‐Seq (Fig EV3 and Table EV1). Around 10,000 genes had an average expression of more than 3 RPM in both D15 and D30 (Fig EV3) and were screened for oscillatory expression over 48 h as determined by JTK‐cycle 50. This revealed 643 and 757 oscillating transcripts (P < 0.05) at D15 and D30, respectively (Fig 4A and Table EV2). The oscillatory transcripts of D15 could result from a starting clock as indicated by small circadian Per2‐dLuc signals at this time point (Fig 3C), but are mostly distinct from the CCGs that were found at D30 (Fig 4B). This limited fraction of overlap might be an underrepresentation, as detecting oscillatory transcription of genes has been shown to rely strongly on sequencing depth 51. Relatively low expression also explains the absence of core clock genes from the rhythmic transcripts (Fig EV5D). Indeed, in our data for both D15 and D30, oscillatory genes had on average more coverage than non‐oscillatory transcripts (Fig EV5E) and shared oscillators between D15 and D30 (n = 80) had higher expression levels than stage‐specific oscillators (Fig EV5F). A fraction of the oscillators (D15 only, D30 only and shared) was found to overlap known mouse cardiac CCGs 14 (Fig 4C) including genes with a known important role in cardiomyocytes (COL4A1, SPON2, SLC23A2, AQP1, and STC1; Fig 4D) 52, 53, 54, 55. These data thus contain common rhythmically expressed clock‐controlled genes between mouse hearts and human ES cell‐derived cardiomyocytes.

Figure 4

Identification of oscillatory transcripts at D15 and D30

643 and 757 oscillators for D15 and D30 cultures, as analyzed using JTK‐cycle (adj. P < 0.05). Heatmaps represent z‐normalized RPM values of the average of three independent replicates. Oscillatory genes were ranked by their phase of expression and visualized using Java TreeView.

Venn diagram of JTK‐cycle detected oscillators for D15 and D30.

Fraction of JTK‐cycle detected oscillators that were previously found to be rhythmically expressed in mouse hearts 14.

Examples of overlapping oscillators between D15, D30 cardiac cells and mouse hearts. Average log2 RPM read counts of three replicates, smoothened over 2 time points ± s.e.m., were plotted. Significance of rhythmicity across 48 h was analyzed using the JTK‐cycle algorithm (* JTK P < 0.05, ** JTK P < 0.005, and *** JTK P < 0.0005).

STRING protein–protein analysis 56 on D30 oscillators that were also found to have circadian expression in mouse hearts (n = 135; Table EV2) revealed a putative, highly interactive network (Fig 5A). D15 rhythmic transcripts that overlap mouse heart oscillators (n = 98), however, did not show such interactions (Fig 5B). Gene ontology (GO) analysis for these oscillators showed enrichment for extracellular matrix formation terms at D15, while D30 oscillators were enriched for cardiac development and stress response terms (Fig 5A and B). Interestingly, the D30 interaction network was centered around UBC (ubiquitin C) (Fig 5A), one of the four genes encoding for ubiquitin in mammals and one of the most abundant proteins in eukaryotic cells 57. Although Ubc is expressed in multiple tissues in mice (http://biogps.org/) 58, it has only been shown to oscillate in the murine heart 52 and (skeletal) muscle 14 (JTK, P = 3.32E‐6 and P = 5.97E‐7, respectively; Fig 5C). This suggests that Ubc is a heart‐ and muscle‐specific CCG in vivo, and concurs with our identification of UBC as a circadian CCG in in vitro D30 human ES cell‐derived cardiomyocytes (D30, JTK, P = 0.0032; Fig 5D). Among the putative UBC interacting partners, several genes were known oscillators in the murine heart according to the CircaDB database (http://circadb.hogeneschlab.org/) 11, 14, 59. Interestingly, many of the oscillating UBC interaction partners in D30 human ES cell‐derived cardiomyocytes were involved in cardiac function (PLN) 60, stress response (BNIP3, RRAGA, DNAJA1 and HSPH1) 61, hypertrophy (RGS2) 62, and even contained therapeutic targets such as TSPO (Translocator protein) 63 (Fig 5D). This indicates that the oscillators that were identified here possibly contribute to multiple molecular mechanisms with a circadian clock dependency, but could also suggest a role for circadian processes in pathophysiological events such as ischemic damage after myocardial infarction.

Figure 5

Circadian stress network results in time‐dependent apoptotic response of human ES cell‐derived cardiac cultures

STRING interaction network of overlapping oscillators between D30 cardiac cells and mouse hearts with corresponding GO analysis. Genes that oscillate at both D15 and D30 are depicted in pink. Mouse heart oscillators were deducted from Zhang et al 14.

Same analysis as in (A) for D15 oscillators.

Ubc mRNA oscillation in mouse hearts and skeletal muscle. Data were obtained from CircaDB (http://circadb.hogeneschlab.org/, deducted from 14, 83). Corresponding CircaDB P‐values are indicated.

Expression levels of four D30 oscillatory genes (UBC, RGS2, PLN, and TSPO) across 48 h. Average log2 RPM CEL‐Seq counts of three replicates were plotted for D15 and D30 and smoothened over two time points ± s.e.m. Significance of rhythmicity across 48 h was analyzed using the JTK‐cycle algorithm (ns: not significant; * JTK P < 0.05, ** JTK P < 0.005, and *** JTK P < 0.0005).

Apoptosis, as measured by caspase 3/7 activity, after doxorubicin administration in D15 and D30 human ES cell‐derived cardiac cells across all samples from (F) and (G). Bottom and top of the boxes are the 25th and 75th percentiles. The line within the boxes represents the median and whiskers denote the minimum and maximum values. The effect of doxorubicin versus DMSO was tested using a Mann–Whitney U‐test (***P < 0.0005).

Apoptosis measured with 6‐h intervals in synchronized D15 cultures after administration of doxorubicin (orange) and DMSO as control (black) for 6 h. Data are represented as the mean ± s.e.m. of three independent replicates. Significance of rhythmicity across 48 h was analyzed using the RAIN algorithm and is indicated (ns: not significant).

Same as in (F) for D30 cultures with the doxorubicin response depicted in red (RAIN, ns: not significant and *P < 0.05). Data are represented as the mean ± s.e.m. of three independent replicates.

Human ES cell‐derived cardiomyocytes show rhythmicity in doxorubicin‐induced apoptosis

Mouse hearts show circadian rhythmicity in their tolerance to ischemia and reperfusion after myocardial infarction 31. In humans, a similar time of the day pattern in the onset and severity of myocardial infarction has been described 25, 32, 33, 34, 64, 65. The combination of time‐dependent pathophysiology and the enrichment of oscillating stress‐associated genes (around UBC) in our CEL‐Seq datasets prompted us to assess whether in vitro derived cardiac cells would show a functional circadian reaction to induced stress. The anthracycline doxorubicin is a widely used anti‐cancer drug that is often administered in the clinic, but is also known to have severe cardiotoxic side effects 66, 67. These effects are recapitulated in in vitro human ES cell‐derived cardiomyocytes in which doxorubicin is known to induce apoptosis and has proven to be a good model for induced cardiotoxicity 68, 69. To determine whether the sensitivity of human ES cell‐derived cardiac cells to doxorubicin‐induced apoptosis displays an oscillatory response, we synchronized cultures at D15 and D30 and administered doxorubicin (10 μM) every 6 h over the course of 48 h (see Materials and Methods). A marked induction of apoptosis was found at both stages, as indicated by elevated active caspase 3/7 levels over control DMSO‐treated samples (Mann–Whitney U‐test, P < 0.0005 for D15 and D30; Fig 5E), with matured D30 cells being more sensitive to doxorubicin than D15 cultures (P < 0.05). Interestingly, the strength of the apoptotic response demonstrated a significant circadian pattern at D30, but not D15 (RAIN, P < 0.05 and P = 0.85 for D30 and D15, respectively; Fig 5F and G), which reveals the functional consequences of a circadian clock in cardiomyocytes. These results highlight the potential of reducing cardiotoxic side effects by the use of time‐based cancer therapy, but also indicate that taking diurnal rhythmicity into account could possibly improve other treatment strategies.

Discussion

Circadian rhythmicity is crucial to heart function, but also influences pathophysiology as indicated by, for instance, diurnal rhythmicity of cardiac damage after infarction 25, 31, 32, 33, 34, 64, 65. As human ES cell‐derived cardiomyocytes are emerging as a powerful tool to model developmental and disease processes as well as being a potential cellular source for regenerative therapies, we examined the presence and possible implications of a functional clock in human ES cells and their cardiac derivatives. While human ES cells do express core clock genes, no circadian clock was observed. Upon differentiation toward cardiomyocytes however, a functional core clock pathway was gradually established (Fig 6) as determined by robust anti‐phasic oscillations of BMAL1 and PER2. This work is the first demonstration of a functional clock in human ES cell‐derived (cardiac) cells and may serve as a paradigm for the emergence of diurnal rhythms in other human pluripotent stem cell‐derived cell types. At D30, 757 CCGs were identified, 18% of which are known to oscillate in the murine heart. Importantly, our data uncover additional transcripts with specific oscillatory behavior in human ES cell‐derived cardiomyocytes. As some of these newly identified oscillators are known to play an important role in human heart physiology (PLN, KCNE4, TSPO, CAV1, RGS2), this stresses the importance of using human cells for modeling cardiovascular processes and disease.

Figure 6

Schematic representation of the emergence and consequences of circadian rhythmicity during directed cardiac differentiation of human ES cells

hESCs: human ES cells. CMs: cardiomyocytes.

Importantly, a defined set of the oscillators could clearly be linked to stress response, which was confirmed by a time‐dependent response to doxorubicin administration. This highlights the possible beneficial effects of drug administration at a specific time of the day to decrease cardiotoxic side effects. Notably, next to explicit clock synchronization steps, such as with forskolin or dexamethasone, simple medium changes can also reset the internal clock of cell cultures 70. Our results demonstrate that circadian mechanisms can influence cellular response to external stressors and thus is an important factor to consider when interpreting experimental results. Our data stress the importance of testing compounds in a time‐controlled manner when using in vitro cultured cardiomyocytes, and may also extent to other ES cell‐based disease models.

Materials and Methods

ESC culture and cardiomyocyte differentiation

NKX2‐5‐eGFP human ES cells 42 (stable reporter line generated from wild‐type HES3 cells 71) were cultured in Essential 8™ medium (Gibco) on Matrigel (BD, Corning) without penicillin/streptomycin. Cells were differentiated in a monolayer toward cardiomyocytes as previously described 45. In short, human ES cells were cultured in E8 until 60% confluent. Cells were then supplemented with 1% DMSO enriched E8 medium for 24 h. On D0, cells were put on BPEL medium supplemented with Activin A (20 ng/ml, R&D Systems), BMP4 (20 ng/ml, R&D Systems), and CHIR99021 (1.5 μM, Axon Medchem). At D3, medium was changed to BPEL with XAV939 (5 μM, Tocris), and on D6, BPEL without any supplements was used. BPEL medium: IMDM, no phenol red (Gibco) and F12 Ham's F12 nutrient Mix (Gibco) in a 1:1 ratio supplemented with 5% (v/v) PFHM‐II (Gibco), 0.25% (w/v) BSA, 1% (v/v) Chemically Defined Lipid Concentrate (Gibco), 0.1% ITS‐X (Gibco), 450 μM α‐MTG (Sigma), 2 mM GlutaMax, 50 μg/ml L‐ascorbic acid 2‐phosphate (Sigma), and 0.25% penicillin/streptomycin (10,000 U/ml, Gibco).

Lentiviral constructs and transduction

Lentiviral plasmids harboring luciferase reporters of the Per2 and Bmal1 promoters were described previously and kindly provided by Prof. Dr. Liu 36, 37, 72. Viral particles were concentrated via ultracentrifugation after three harvests in HEK293T cells. Cells were transduced with concentrated Bmal1‐dLuc or Per2‐dLuc lentivirus 2 days before circadian bioluminescent measurements.

Bioluminescent recording and data analysis

Human ES cell‐derived cardiac cells were differentiated for up to 45 days and transduced with lentiviral reporters at described time points after synchronization with 100 nM dexamethasone 47 for 2 h. Subsequently, medium was changed to recording medium [BPEL, 10 mM HEPES, 100 μM D‐Luciferin Potassium Salt (Promega)]. Human ES cells were cultured in E8 medium and synchronized for 2 h with forskolin 39. Forskolin was chosen as a synchronizing agent for human ES cells, since dexamethasone has been implemented in multiple stem cell differentiation protocols and might therefore induce premature differentiation 73, 74, 75, 76, 77, 78. Subsequently, medium was changed to ES recording medium (E8, 10 mM HEPES, 100 μM D‐Luciferin Potassium Salt (Promega)). Culture dishes were sealed with high vacuum grease (Dow Corning) and monitored via the use of a LumiCycle32 device (Actimetrics) at 37°C. Bioluminescence from each dish was continuously recorded (integrated signal of 70 s with intervals of 10 min). Raw data (counts/s) were baseline subtracted (polynomial order 3).

Microscopic real‐time bioluminescence analysis

Human ES cells were differentiated for up to 30 days and transduced with Bmal1‐or Per2‐dLuc lentivirus 2 days before recording. Bioluminescence was assessed with an LV200 microscope (Olympus) in a humidified chamber under 5% CO2, at 37°C. Bioluminescence was detected for multiple consecutive days, using an EM CCD camera (Hamamatsu), with exposure times of 1 h. Image series were analyzed in ImageJ. Cells were synchronized with 100 nM dexamethasone for 2 h and changed to normal BPEL medium, containing 100 μM D‐Luciferin Potassium Salt (Promega). For pure cardiomyocyte population experiments, human ES cell‐derived cardiac populations were sorted with a FACS Jazz flow cytometer (BD Biosciences) based on GFP positivity, replated on Matrigel‐coated dishes, and bioluminescence was assessed 7 days later.

Immunostaining

Cells were fixed with 4% paraformaldehyde (PFA) for 15 min, blocked for 1 h in blocking buffer (5% FBS, 0.25% Triton X‐100 in PBS), and stained for OCT4 (SantaCruz, #5279), NANOG (Cell Signaling, #3580S), TNNT2 (ThermoFisher Scientific, #MA5‐12960), ACTN2 (Sigma, #A7811), MEF2C (Cell Signaling, #sc‐13266), GFP (Abcam, #ab6556) in staining buffer (1% BSA, 0.25% Triton X‐100 in PBS). Nuclei were stained with Hoechst for 15 min. Images were made using a spinning disk microscope (PerkinElmer).

RNA isolation and CEL‐Seq

Cardiac cells were derived from human ES cells in 48‐well plates. After synchronization, biological triplicates (independent wells) with comparable cardiac purity were collected every 4 h over the course of 48 h (ZT4‐ZT48) (Fig EV3). RNA was extracted using the standard TRIzol (Invitrogen) protocol, and 10 ng of total RNA per sample was used for library preparation and sequencing. RNA was processed as described previously 49, 79, and paired‐end sequencing was performed on the Illumina Nextseq platform with a read length of 75 base pairs. Read 1 was used to identify the sample barcode and library index, while read 2 was aligned to the hg19 human RefSeq transcriptome (downloaded from the UCSC genome browser) using BWA 80. CEL‐Seq only sequences the most 3′ end of a transcript, generating one read per transcript. Reads that mapped equally well to multiple locations were discarded. Around 500,000 reads were sequenced per sample. Samples were reads per million (RPM) normalized (Table EV1).

Quantitative RT–PCR

Purified RNA was treated with DNAse (Promega) and reversibly transcribed with Superscript III reverse transcriptase (ThermoFisher Scientific). qRT–PCR on biological triplicate samples was carried out in triplicate (technical replicates) in CFX‐384 Touch™ Real‐time PCR detection system (Bio‐Rad). PPIA was used as housekeeping gene, and fold changes were calculated to the lowest values among all replicates. Primer sequences: PPIA (fw): ttctgctgtctttgggacct, PPIA (rv): caccgtgttcttcgacattg, NANOG (fw): cagccctgattcttc, NANOG (rv): tgcatctgctggaggctgag, POU5F1 (fw): ctgaagcagaagaggat, POU5F1 (rv): gggccgcagcttacacat, ISL1 (fw): ctgcttttcagcaactggtca, ISL1 (rv): ggactggctaccatgctgtt, NKX2‐5 (fw): caagtgtgcgtctgcctttc, NKX2‐5 (rv): ctttcttttcggctctagggtcct, ACTC1 (fw): atgccatcatgcgtctggat, ACTC1 (rv): acgttcagcagtggtgacaa, KCNJ2 (fw): tgggtcttgggaattctggttt, KCNJ2 (rv): gaacatgtcctgttgctggcg, SERCA2A (fw): cgaacccttgccactcatct, SERCA2A (rv): ccagtattgcaggttccaggt, BMAL1 (fw): ggctcatagatgcaaaaactgg, BMAL1 (rev): ctccagaacataatcgagatgg, PER2 (fw): ggccatccacaaaaagatcctgc, PER2 (rv): gaaaccgaatgggagaatagtcg, CRY1 (fw): ctccatgggcactggtctcagtg, CRY1 (rv): tccccaccaatttcagctgcaac, CRY2 (fw): ccaagagggaagggcagggtagag, CRY2 (rv): aggatttgaggcactgttccgagg, CLOCK (fw): aagttagggctgaaagacgacg, CLOCK (rv): gaactccgagaagaggcagaag, NR1D1 (fw): acagctgacaccacccagatc, NR1D1 (rv): catgggcataggtgaagatttct.

Western blotting

Cells were lysed in RIPA buffer, and protein concentration was measured using a BCA assay (ThermoFisher Scientific). 12.5 μg protein lysate was loaded, separated by 10% SDS–PAGE, and transferred to a nitrocellulose membrane. Membranes were blocked with 5% milk powder (Nestlé) in T‐BST and probed with anti‐BMAL1 (1:1,000, #ab3350, Abcam), anti‐CRY1 (1:1,000, #13474‐1‐A, Proteintech), or anti‐CLOCK (1:250, #PA1‐520, ThermoFisher Scientific) antibodies, followed by a peroxidase‐conjugated antibody (1:5,000, #sc‐2004, Santa Cruz). ECL Plus Western blotting substrate (#32132, ThermoFisher Scientific) was used for chemiluminescence detection with an ImageQuant™ LAS 4000 imager (GE Healthcare). HRP‐coupled anti‐β‐ACTIN (1:5,000, #5125S, Cell Signaling) was used as a loading control. Band intensities were calculated with ImageJ.

Apoptosis measurements

Human ES cells were differentiated in 96‐well white walled plates for the course of 15 and 30 days. Cardiac cells were synchronized with 100 nM of dexamethasone for 2 h, and 10 μM of doxorubicin HCl (Sigma D1515) was administered at 6‐h intervals for a total time of 6 h. Apoptosis levels of three replicate wells (per condition), represented by active caspase 3 and caspase 7 levels, were measured using a CaspaseGlo 3/7 kit (Promega) following manufacturer's instructions. Bioluminescence was read out with a Centro microplate luminometer (Berthold Technologies).

JTK‐cycle analysis

RPM‐normalized read counts were obtained for each sample. As lowly expressed genes are typically not picked up robustly using CEL‐Seq, genes with an average of > 3 RPM across all time points (ZT4‐ZT48) as well as the replicates were selected for JTK‐cycle analysis. Around 10,000 genes reached this threshold in both D15 and D30, and these form the list on which JTK‐cycle was run (Table EV2). The following settings were used: jtkdist (12,3), periods (6:6), jtk.init (periods, 4). Significant oscillators with an adjusted P‐value of < 0.05 were selected for further analyses. To identify mouse heart oscillators, JTK was run with similar settings on normalized GC‐RMA intensity values of 24 samples (CT18‐CT62, sampled every 2 h) for 35,556 genes downloaded from the GEO‐database 81 (accession GSE54652) 14.

STRING and gene ontology analysis

STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) database (www.string-db.org) was used to investigate the relationship between the overlap of known murine heart oscillators and identified D15 and D30 oscillators 56. Gene ontology terms were retrieved via www.string-db.org.

CircaDB gene expression website

The circadian expression database (CircaDB, http://circadb.hogeneschlab.org/) is an open access online platform 59 compiling circadian gene expression profiles from microarrays and RNA sequencing experiments 11, 14, 18, 82, 83, 84, 85, 86, 87, 88. The embedded JTK‐cycle algorithm defines the significance of rhythmic gene expression.

Statistics

All data were shown as means ± s.e.m. Student's t‐tests were carried out to assess differences between qRT–PCR mean values within the same experiments. One‐way ANOVA, followed by a Bonferroni post hoc test, was carried out to test increasing mRNA levels of clock factors during directed cardiac differentiation. A difference of P < 0.05 was considered significant. To calculate general induction of apoptosis upon doxorubicin, a Mann–Whitney U‐test was used. Differences in doxorubicin‐effect sizes between D15 and D30 cardiac cells were assessed via non‐overlapping 95% effect interval sizes. Statistical analyses to detect circadian oscillations in RNA levels (qRT–PCR) as well as doxorubicin‐based apoptosis were performed by RAIN 40.

Data availability

Primary data

Dierickx P, Vermunt MW, Muraro MJ, Creyghton MP, Doevendans PA, van Oudenaarden A, Geijsen N, Van Laake LW (2017) Circadian networks in human embryonic stem cell‐derived cardiomyocytes. Gene Expression Omnibus GSE97142.

Referenced data

Zhang R, Lahens NF, Ballance HI, Hughes ME, Hogenesch JB (2014) A circadian gene expression atlas in mammals: Implications for biology and medicine. Gene Expression Omnibus GSE54652.

Author contributions

PD, PAD, LWVL, and NG conceived and designed the research; PD conducted and analyzed the experiments; MJM supervised by AvO, made sequencing libraries, and processed CEL‐Seq data; PD and MWV supervised by MPC, analyzed and interpreted the CEL‐Seq data; PD, MWV, NG, and LWVL wrote the paper.

Conflict of interest

The authors declare that they have no conflict of interest.

Expanded View Figures PDF

Click here for additional data file.

Table EV1

Click here for additional data file.

Table EV2

Click here for additional data file.

Review Process File

Click here for additional data file.