Differentiated derivatives of human pluripotent stem cells (hPSCs) are often considered immature because they resemble foetal cells more than adult, with hPSC-derived cardiomyocytes (hPSC-CMs) being no exception. Many functional features of these cardiomyocytes, such as their cell morphology, electrophysiological characteristics, sarcomere organization and contraction force, are underdeveloped compared with adult cardiomyocytes. However, relatively little is known about how their gene expression profiles compare with the human foetal heart, in part because of the paucity of data on the human foetal heart at different stages of development. Here, we collected samples of matched ventricles and atria from human foetuses during the first and second trimester of development. This presented a rare opportunity to perform gene expression analysis on the individual chambers of the heart at various stages of development, allowing us to identify not only genes involved in the formation of the heart, but also specific genes upregulated in each of the four chambers and at different stages of development. The data showed that hPSC-CMs had a gene expression profile similar to first trimester foetal heart, but after culture in conditions shown previously to induce maturation, they cluster closer to the second trimester foetal heart samples. In summary, we demonstrate how the gene expression profiles of human foetal heart samples can be used for benchmarking hPSC-CMs and also contribute to determining their equivalent stage of development.
In the developing embryo, the heart is one of the first organs to be fully formed. It develops from a linear tube into a four-chambered organ through a complex looping process that leads to the formation of the ventricles, the atria and the outflow tract. In humans, the heart starts beating around 6 weeks of gestation and pumps blood though a closed circulatory system to provide nutrients and oxygen and remove waste products from organs as they develop. Later, the septum separates the atria into left and right halves and closes completely after birth. The left and right ventricles are separated before birth and their walls develop into strong muscles (Moorman et al., 2003). The left ventricular wall is thicker than the right because it pumps oxygenated blood from the lungs to all parts of the body via the aorta. The right side of the heart receives de-oxygenated blood and pumps it through the lungs to re-oxygenate. Although much is known about the molecular mechanisms that drive heart formation and morphogenesis in laboratory animals (Harvey, 2002), little equivalent data is available on the human heart. This is important for understanding how specific mutations in different genes (i.e. missense mutations), rather than knockouts commonly used in experimental animals, affect human heart development and function as well as validating models of hereditary heart disease based on patient-derived human induced pluripotent stem cells (hiPSCs).
The four chambers of the mammalian heart express different genes that determine their physiological properties (Small and Krieg, 2004). Most studies to date have analysed the transcriptome of adult human atria and ventricles (Asp et al., 2012), or have performed transcriptional analysis on human auricle (part of the atrium) removed during normal clinical procedures on diseased adult hearts (Sanoudou et al., 2005; Nanni et al., 2006). Collection of healthy human ventricular tissue is more difficult and is usually only available when donor hearts are not used for transplantation. Here, we collected tissue from the chambers of human foetal hearts in isogenically matched combinations of atria and ventricles from either the first (T1) or second (T2) trimester of development and examined gene expression as a function of gestational age and region of the heart. Using the advanced microarray techniques now available, we were able to compare transcriptomes of foetal tissue using small amounts of RNA (<0.5 μg). These studies not only provided better insight into how the human heart develops, but also allowed us to compare these samples with cardiomyocytes (CMs) derived from human pluripotent stem cells (hPSCs).
hPSCs are now frequently used to investigate CM differentiation in the early embryo and also to model inherited cardiac disease (Davis et al., 2011). Although it is widely acknowledged that hPSC-derived CMs (hPSC-CMs) are immature (Yang et al., 2014), how they compare with the individual chambers of human foetal hearts is not clear. We have examined hPSC-CMs cultured in conventional (defined) differentiation medium (van den Berg et al., 2015), and in medium recently reported to induce functional maturation as evidenced by increased force of contraction (Ribeiro et al., 2015). We found that the foetal gene expression profiles allowed us to estimate the maturity of the hPSC-CMs.
RESULTS AND DISCUSSION
Global gene expression analysis of foetal and adult heart samples show changes with age
We determined the transcriptional profiles of human foetal hearts collected during T1 and T2 of normal human development. We separated T1 hearts into atrium and ventricle. As T2 hearts were larger, we collected the four chambers separately. We also included a commercially available sample of pooled adult hearts as a common reference for normalization of future samples (Fig. 1A; supplementary material Table S1). Human foetal heart samples were a mixed population of cardiomyocytes, fibroblasts and endothelial cells because there are no specific cell surface antibodies for cardiomyocytes that would allow them to be sorted from primary heart tissue. However, the majority of the cells in the foetal heart are cardiomyocytes, which decreases postnatally when cardiomyocyte division ceases (Bergmann et al., 2015).
To compare the variability between all samples, we performed gene expression cluster analysis (Fig. 1B). We found that the foetal heart samples showed distinct gene expression profiles between T1 and T2. From the T1 and T2 samples, the atria and ventricles clustered separately, except for the foetal heart (FH) FH5-right atrium (RA), FH5-right ventricle (RV) and FH5-left ventricle (LV) samples, which clustered together as individual number 5. T2 ventricle samples clustered closer to the adult heart reference samples, possibly because the contribution of ventricles to the pooled adult reference sample is greater than the (smaller) atria.
We also grouped samples based on age (7 weeks, 15 weeks, 20 weeks and adult) and investigated all differentially expressed genes (DEGs) in the first and second trimester and adult hearts. We identified eight distinct transcriptome clusters (Fig. 1C). The DEGs in Cluster 1 were upregulated at 7 weeks, showed no changes in 15- and 20-week samples and decreased in adult heart. Cluster 2 showed a similar pattern; DEGs in both clusters were involved in cell cycle regulation and chromatin organization. DEGs in Cluster 3 were upregulated during foetal development and downregulated in adult heart. This cluster included genes important in development, cell division and matrix organization that are less important in adult hearts unless the heart has been damaged, for example by myocardial infarction. As the heart develops and ages, CM proliferation, which is essential during early heart development, decreases (Bergmann et al., 2015). Cluster 4 showed genes that gradually decreased over time and were involved in cell cycle regulation and chromatin organization. Cluster 5 contained DEGs involved in gaseous substance transport; these were downregulated at week 7 and in adult heart but were upregulated at intermediate stages. Clusters 6, 7 and 8 consisted of DEGs that increased over time and are important in metabolic processes, muscle organization and contraction. The foetal heart depends on carbohydrate synthesis, but it also prepares for the switch in metabolism to fatty acid oxidation shortly after birth (Taegtmeyer et al., 2010). Gene ontology terms for biological processes (GO.BP) for each cluster are listed in supplementary material Table S2.
First and second trimester atria and ventricles have distinct gene expression signatures
To investigate genes that are important for development of atria and ventricles, we divided samples into four groups according to their origin and age: atria T1 (A1), ventricles T1 (V1), atria T2 (A2) and ventricles T2 (V2). Four comparisons were made to investigate differences in gene expression between atria and ventricles, and also between T1 and T2. Gender did not appear to influence gene expression in the heart with only genes located on the sex chromosomes (i.e. XIST, RPS4Y2, DDX3Y, RPS4Y1, EIF1AY) being differentially expressed between age-matched male and female samples. In subsequent analyses, we therefore combined male and female samples. Using an absolute log2 fold difference ≥1.5 in combination with a significance threshold of P-adjusted value of 0.05, we identified a total of 156 DEGs. Two-way cluster analysis of all DEGs revealed distinct transcription profiles within all four groups (Fig. 2A). In contrast to the clustering based on the global gene expression in Fig. 1B, all samples now separated based on the trimester and chamber subtype. We found 24 DEGs in T1 atrium versus ventricle (A1/V1), 34 in T2 atrium versus ventricle (A2/V2), 110 in T2 versus T1 atria (A2/A1) and 39 in T2 versus T1 ventricles (V2/V1) (supplementary material Table S3). Fig. 2B shows the overlap between the four comparisons, with the upper Venn diagram focusing on differences between atria and ventricles and the lower diagram on age. Overall, few genes overlapped between atria and ventricles and DEGs included genes that were reported to be expressed in a chamber-specific pattern in both mouse and human, such as NR2F1, MYL2 and KCNA5. Between T1 and T2, the number of overlapping genes per chamber was higher and correlated to chromatin and nucleosome structure (downregulated) and extracellular matrix and collagen organization (upregulated). The volcano plots in Fig. 2C display the DEGs, with selected genes highlighted based on the results here and from earlier publications on the adult cardiac transcriptome (Ng et al., 2010; Lu et al., 2014).
In both T1 and T2 atria versus ventricle, MYL2 was downregulated, reflecting its importance in ventricle contraction, whereas KCNA5, encoding the potassium channel Kv1.5, which conducts the ultra-rapid activating delayed rectifier K+ current (IKur), a major repolarizing current in human atria (Christophersen et al., 2013), was upregulated. MYL3 and MYH7 were also described previously as differentially expressed in gene and protein studies comparing atria and ventricles in adult (Asp et al., 2012; Lu et al., 2014). Among the genes expressed at higher levels in T1 and T2 atria compared with similarly aged ventricles were NR2F1 (also known as COUP-TFI) and RELN. NR2F1 was recently shown to be enriched in the atria of human foetal as well as adult hearts and to regulate atrial-specific ion channel genes in atrial hPSC-CMs (Devalla et al., 2015). RELN is involved in neuron migration and brain development (Tissir and Goffinet, 2003), but has not previously been described in foetal heart development. Other studies have also detected higher expression levels of RELN in human adult atria (Kaab et al., 2004), although its exact function in the heart is unknown.
Chromatin remodelling and histone modifications are also known to have an important role in heart development (Han et al., 2011; Chang and Bruneau, 2012). Genes encoding histones that influence nucleosome structure and are important in compaction of DNA (HIST1H3I, HIST1H2BM, HIST1H2AI) were mainly downregulated in T2 ventricles and atria, and to our knowledge have not previously been described in cardiac development. Genes important for extracellular matrix and collagen fibril organization (COL1A2, COL12A1, COL15A1, DPT) were upregulated in the T2 samples, indicating further development of the cardiac scaffold and maturation of the heart. Among the DEGs were also genes involved in cardiac development, electrical currents and sarcomere structure, such as BMP10, APLNR, PLN, TNNI3K and MYOM2. BMP10 has been reported in mouse heart development (Chen et al., 2004) and previous studies have also detected BMP10, which is involved in the trabeculation of the heart, to be more highly expressed in atria (Asp et al., 2012; Lu et al., 2014).
Selected gene ontology terms show chamber-biased expression
To explore the functional characteristics of the DEGs, we performed GO.BP. Fig. 3 and supplementary material Table S4 display GO.BP terms with adjusted P-value ≤0.01. As expected, among the genes expressed at lower levels in A2 versus V2 were those related to ventricle development, contraction and muscle morphogenesis and structure (Fig. 3A). In A2 compared with A1, chromatin and nucleosome organizational genes were downregulated, suggesting that a process of active chromatin remodelling is slowing down at this stage of development. No GO.BP terms were enriched when V2 and V1 were compared. GO.BP terms involved in extracellular matrix organization, wound healing and blood coagulation were upregulated in T2 compared with T1 and included genes such as VWF and APLNR (Fig. 3B). Foetal genes are typically upregulated during remodelling of the adult heart, for example after myocardial infarction, owing to the activation of pathways such as wound healing and stress responses (Gidh-Jain et al., 1998; Zgheib et al., 2014). We also found a significant over-representation of genes in atrial samples involved in neuron generation, forebrain development or neuron migration, such as NR2F1 and NR2F2. This is likely to be due to innervation of the heart and the control of cardiac rhythm by the autonomous nervous system (the vagus nerve around the sinus node) at this stage of development (Kimura et al., 2012), or due to genes that are expressed both in atria and the brain (Alfano et al., 2014).
The foetal heart transcriptome can indicate the maturation state of hPSC-CMs
To demonstrate the utility of this dataset, we compared human embryonic stem cell (hESC)- and hiPSC-derived CMs with the foetal heart samples with the expectation that this would provide insight into how their gene expression profiles relate to primary cardiac tissue (Fig. 4A). An earlier study described the similarities and differences between foetal heart samples and CMs but the hESC-derived CMs included were derived from mixed population clusters and foetal heart samples were only obtained from third trimester donors (Synnergren et al., 2012). Here, we had the opportunity to compare the hPSC-CMs with primary tissue at earlier stages of foetal development and from different chambers of the heart.
hPSC-CMs were initially purified from mixed populations of differentiated cells maintained in LI-BPEL medium on the basis of the reporter gene eGFP, which had been inserted into the cardiac transcription factor NKX2-5 (NKX2.5-GFP) genomic locus of hESCs (Elliott et al., 2011) and hiPSCs (C.W.v.d.B., C.L.M. and R.P.D., unpublished). The differentiation protocols used here yield primarily ventricular-like cardiomyocytes on the basis of action potentials in patch clamp electrophysiology (Davis et al., 2012; Bellin et al., 2013) and by excluding the NKX2.5-GFP− cardiomyocytes, we also excluded pacemaker-like cells (Birket et al. 2015). Furthermore, we also examined hESC-derived CMs that were cultured in commercially available maturation medium (MM) containing T3 hormone as a principal component (Ribeiro et al., 2015). To examine similarities and differences between in vitro hPSC-CMs and foetal heart, we compared transcriptional profiles of the samples with the foetal heart from each trimester. Initially, we performed hierarchical cluster analysis of all PSC-derived cardiomyocytes together with all foetal heart samples (supplementary material Fig. S1). As expected, all the foetal heart samples clustered closer to each other than to the hPSC-CMs. However, we did observe that hPSC-CMs cultured in MM clustered closer to the foetal heart samples than to the hPSC-CMs maintained in standard culture medium, suggesting that the cardiomyocytes had indeed developed further. To investigate which foetal heart age group the hPSC-CMs in LI-BPEL or MM most closely matched, we performed cluster analysis of the hPSC-CMs with the foetal heart samples from T1 and T2 separately (Fig. 4B,C). The cells maintained in the regular LI-BPEL culture medium were more closely related to T1 foetal heart samples than were hPSC-CMs cultured in MM. When comparing the in vitro CMs to T2 samples, the CMs that had been maintained in MM more closely resembled the T2 heart. The partial maturation of the hESC-CMs in MM that we observed correlated well with another recent study that also showed that individual hPSC-CMs cultured in MM developed functional features closer to that of T2 foetal cardiomyocytes, with most strikingly a greater than twofold increase in contraction stress compared with hPSC-CMs cultured in regular differentiation conditions (Ribeiro et al., 2015).
We have analysed here gene expression profiles in a rare and complete set of isogenic foetal heart samples and described differences in genes expressed between different chambers of the heart during the first and second trimester of development. We showed that microarray analysis could be performed on RNA samples as small as 50 ng, a technical advance that allowed inclusion of hearts at the very earliest stages of development. Our results revealed a group of nucleosome- and histone-related genes expressed in human foetal hearts that to our knowledge have not been described before in cardiac development in mice. Furthermore, we demonstrated how the foetal heart dataset can be used to benchmark hPSC-CMs in terms of their maturation state. The question of how mature hPSC-CMs are frequently arises and our dataset provides a means of answering this by global gene expression rather than on the basis of specific markers. We have included commercially available reference RNA sets in the analysis in order to characterize future sets of hPSC-CMs cultured in conditions that further induce maturation improvements. These reference sets can be used by other laboratories, as well as with other microarray platforms, for normalization and benchmarking of future datasets to the human foetal heart dataset. The foetal heart dataset is provided here as a resource to the community.
MATERIALS AND METHODS
Foetal heart sample collection and ethics statement
Human foetal heart samples were collected from eight healthy individuals after elective abortions at various gestational ages (7, 10, 15, 20 and 20+ weeks of gestation) determined using obstetric ultrasonography based on crown-rump length measurements. The Medical Ethical Committee of the Leiden University Medical Center (protocol 08.087) approved the use of human foetal material and informed written consent was obtained in accordance with the World Medical Association Declaration of Helsinki guidelines.
Total RNA was extracted from the samples using standard isolation techniques. The quality and integrity of the RNA samples was confirmed using Lab-on-Chip RNA 6000 Nano and RNA 6000 Pico (both Agilent) on the Agilent 2100 Bioanalyzer (Agilent Technologies) by ServiceXS B.V. (Leiden, The Netherlands).
Two reference samples (Universal Human Reference RNA, Cat #740000, Stratagene and Human Normal Heart Donor Pool Cat #R1234122-P, lot #A509251, Biochain) were co-hybridized with the experimental samples (supplementary material Table S1). These references were used for normalization between different arrays. For whole-genome microarray of foetal heart samples, biotinylated ss-cDNA was prepared using the NuGEN Ovation PicoSL WTA v2 System (NuGEN) according to the manufacturer's protocol using ∼50 ng total RNA. For hESC- and hiPSC-derived CMs, biotinylated cRNA was prepared using the Illumina TotalPrep RNA Amplification Kit (Ambion) according to the manufacturer's specifications using ∼200 ng total RNA. Hybridization and processing of all samples was performed on Illumina HumanHT-12 v4 microarray chips by ServiceXS B.V. (Leiden, The Netherlands). Microarray data have been deposited in Gene Expression Omnibus under accession number GSE71148.
Data analysis and statistics
The raw intensity values of the microarray data were normalized by variance stabilizing normalization using the vsn R package (Huber et al., 2002) and subsequently normalized by the common references. When a gene had more than one microarray probe, the one with the highest variance across the samples was used for subsequent analysis. The differential expression analysis was performed using the limma R package (Smyth, 2004). Genes were binned into 30 bins by the intensity and the t-test applied to each bin. The P-value was corrected for multiple testing by the Benjamini-Hochberg method with an adjusted P-value cut-off of 0.05. Genes with the mean absolute log2 fold change <1.5 were discarded. Unsupervised hierarchical clustering of all detected genes was performed using the Euclidean distance and complete linkage method. The same parameters were used for clustering the DEGs. Fuzzy clustering of DEGs among foetal and adult heart samples was performed using the R package Mfuzz. The optimal number of clusters was detected using the C-means clustering algorithm implemented in the R package e1071. GO.BP was downloaded from http://www.geneontology.org. Fisher's exact test was performed to identify the statistical enrichment of these categories using the differentially up- or downregulated genes as the test set. All detected genes were taken as the background set. The P-value was corrected for multiple testing by the Benjamini-Hochberg method. Categories with an adjusted P-value <0.01 and odds ratio >1.0 were considered significantly enriched.
Cardiac differentiation of human ESCs and iPSCs
hPSCs were differentiated to CMs as previously described (Elliott et al., 2011; van den Berg et al., 2015) and maintained either in LI-BPEL (Elliott et al., 2011), or according to the manufacturer's protocol (Pluricyte Medium, Pluriomics). To isolate CMs from the NKX2.5-GFP reporter hESC and hiPSC lines (Elliott et al., 2011) (C.W.v.d.B., C.L.M., R.P.D., unpublished), the cells were dissociated (van den Berg et al., 2015) and sorted based on GFP expression using a BD FACSAria III Cell Sorter (BD Biosciences).
The authors thank C. Grandela for experiments and providing hPSC-CM samples; H. D. Devalla for input in early stages of the study; and M. S. Roost for advice on array analysis.
C.W.v.d.B. performed the experiments; C.W.v.d.B. and S.O. performed data analysis; C.W.v.d.B., S.O., A.d.S., R.P., R.P.D. and C.L.M. developed concepts, interpreted the results and prepared the manuscript; S.M.C.d.S.L., L.v.I., S.R.B. and L.G.T. performed experiments and contributed reagents.
This work was funded by the Netherlands Institute of Regenerative Medicine (NIRM) [FES0908 to S.R.B. and C.L.M.]; an AFR Postdoctoral Grant from the Fonds National de la Recherche Luxembourg (FNR) [7682104/PDR to S.O.]; the European Research Council [ERCAdG 323182 STEMCARDIOVASC to C.L.M.]; the Netherlands Heart Foundation (NHS) [CVON-HUSTCARE to C.L.M. and R.P.]; the Netherlands Organisation for Scientific Research (NWO-Aspasia) [015.007.037 to S.M.C.d.S.L.]; and Interuniversity Attraction Poles (IAP) [P7/07 to S.M.C.d.S.L. and C.L.M.].
S.R.B., R.P. and C.L.M. are co-founders and SRB CSO of Pluriomics B.V.