Three hESC lines [H1 (male), 429 (female) and 360 (male)] and four hiPSC lines (established from two independent Detroit 551 (female) clones (clones 7 and 10) and one CRL 2097 (male) iPSC clone, CRL-8) were employed for this study. The details of reprogramming and characterization of the hiPSC lines are published elsewhere (Siller et al., 2015, 2016; Balafkan et al., 2020; Liang et al., 2020). The hPSCs were cultured using standard procedures in a 5% CO₂ incubator at 37°C. hPSCs were maintained in feeder-free conditions in Essential 8 Medium (E8) (Gibco™). The passage number of the hiPSCs was between 23 and 55 for all experiments, and there was no preference for using a specific passage number for an experiment.

Cardiomyocyte differentiation was performed in 96-well microplates as previously described (Balafkan et al., 2020). Briefly, hiPSCs were seeded and propagated on Geltrex (Gibco™), under feeder-free conditions in E8. Within 3 days, when cells reached the optimum confluency (60–70%), cardiomyocyte differentiation was started by applying the GSK3 inhibitor CHIR99021 in Roswell Park Memorial Institute (RPMI) 1640 medium supplemented with B27 without insulin (RPMI-B27), in a cell-concentration-dependent manner. After 24 h, medium was changed to RPMI-B27 without CHIR99021. The differentiation process was continued by adding 5 μM of inhibitor of WNT production-2, IWP2, diluted in RPMI-B27, 72 h post-differentiation induction for 48 h. Fresh RPMI-B27 medium was provided on day 5, and from day 7, cells were fed with fresh RPMI medium supplemented with B27 with insulin without extra supplement, every 2 days. Differentiated cardiomyocytes start beating around day 10, and it becomes more synchronized after day 12 (see Supplementary Movies 1, 2).

MagMAX^TM-96 Total RNA Isolation Kit was used for RNA isolation from cultured cells. EXPRESS One-Step Superscript qRT-PCR Kit (Invitrogen™) was used for cDNA synthesis and real-time PCR using TaqMan probes (Supplementary Table 1) on an Applied Biosystems 7,500-Fast real-time PCR System. All real-time PCR were performed in triplicate, and the average of Ct values was normalized to the geometric mean of ACTB and GAPDH as endogenous control genes, and the result (dCt) used for further analysis. For a detailed description, see the Supplementary Material.

RNA sequencing (RNA-seq) analyses were performed independently for two datasets. The first dataset (Dataset A, N = 22) was composed of either hESC or hiPSC lines collected from S1 (undifferentiated cells at day 0) and S4 (Isl1⁺ progenitor cells at day 7), from three independent differentiation runs of two hESC lines (429 and 360), two independent differentiation runs of Detroit-7 and CRL-8, and one of Detroit-10. The second dataset (Dataset B, N = 12) was composed of samples collected at four stages spanning S1–S5 (specifically days 0, 2, 5, and 15) from H1. RNA-seq sample information is provided in Supplementary Tables 3, 4. Differential gene expression analyses were performed using the DESeq2 R package version 1.26 (Love et al., 2014) with default parameters, using cell line as a covariate in the model for the first dataset. A corrected p-value of 0.05 and log₂ (fold change) of 1 were set as the threshold for significantly differential gene expression. Genes then were tested for enrichment using the gene score resampling method implemented in the ermineR package version 1.0.1, an R wrapper package for ermineJ (Gillis et al., 2010) with the complete Gene Ontology (GO) database annotation (Ashburner et al., 2000), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000) to obtain lists of upregulated and downregulated pathways for each cohort. The source code for the RNA-seq analyses is available in the GitLab repository¹ under GPL public license v3.0. For a detailed description of the methods, refer to the Supplementary Material.

Cells were dissociated into single-cell suspension using TrypLE^TM Express Enzyme (Gibco™) and stained with Zombie Red^TM Fixable Viability Kit (BioLegend, San Diego, CA, United States) according to the manufacturer’s instructions. Single cells were fixed, blocked, and stained accordingly (for a detailed description, see the Supplementary Material). Antibodies are listed in Supplementary Table 1. At least 30,000 events were collected for the target marker using a Sony cell sorter SH800 (Sony Biotechnology Inc., San Jose, CA, United States), and collected data were analyzed and presented by FlowJo V.10.5.0 (FlowJo LLC, Ashland, OR, United States). To conduct inter- and intra-sample analyses from different runs and minimize the batch-to-batch variation, the flow cytometer was calibrated prior to quantitative fluorescence intensity measurements using Quantum^TM Alexa Fluor^® 488 MESF (molecules of equivalent soluble fluorophore; Bangs Laboratories, Fishers, IN, United States), and collected median fluorescence intensity (MFI) was normalized to MESF as an external control. All gates were adjusted according to fluorescence minus one control (FMO).

Cells were seeded either on Geltrex-coated cover slips or in Millicell^® EZ SLIDES (Merck Millipore, Billerica, MA, United States), fixed, and stained accordingly (for a detailed description, see the Supplementary Material). A final concentration of 10 μg/ml was used for all primary antibodies, and they were incubated overnight at 4°C. Confocal microscopy images were taken on either a Zeiss LSM 510 META (Carl Zeiss, Oberkochen, Germany) or a Leica TCS SP5 (Leica Microsystems GmbH, Wetzlar, Germany) at the Molecular Imaging Center (MIC), University of Bergen; and data analysis and image editing were done with Fiji (Schindelin et al., 2012). Antibodies are listed in Supplementary Table 1.

Samples were prepared by the MIC, Department of Biomedicine, University of Bergen; and grids were imaged using a Jeol JEM-1230 transmission electron microscope (TEM) at 80 kV. For a detailed description, see the Supplementary Material.

MtDNA quantification and deletion assessment were performed on DNA isolated from cultured cells using MagMAX^TM-96 DNA Multi-Sample Kit (Invitrogen™) and real-time PCR, as well as long-range PCR, as previously described (Tzoulis et al., 2013). A commonly deleted region (MT-ND4) in the major arc of mtDNA and a rarely deleted region (MT-ND1) were utilized to quantify deletion, and MT-ND1 was compared with amplification of a single-copy nuclear gene (APP) to assess the number of mtDNA copies. For a detailed description (see Supplementary Material).

Respiration and acidification rates were measured on a monolayer culture of undifferentiated hPSCs and cells at S5 using a Seahorse XF^e96 extracellular flux analyzer (Agilent, Santa Clara, CA, United States). For a detailed description (see Supplementary Material). To correct the final results for differences in cell size and number between the undifferentiated hPSCs and differentiated cells at S5, we measured the total protein concentration for each well using absorbance at 280 nm. All results were reported as pmol O₂ per min after normalization to total protein concentration. The XF reader software (Wave Desktop 2.4) was used to analyze the data.

We used a previously described protocol (Rowe and Boletta, 2013) to analyze the mitochondrial membrane potential (ψm) independent of cell volume using tetramethylrhodamine methyl ester (TMRM). We quantified the MFI of TMRM before and after applying the uncoupler carbonyl cyanide 4-(trifluoromethoxy)phenylhydrazone (FCCP) to dissipate the mitochondrial membrane potential. The difference between geometric means for TMRM- and TMRM–FCCP-treated cells provided the membrane potential normalized to cell size. In addition, to normalization for mitochondrial mass differences, we calculated the ratio of median MFI of TMRM to MFI of TOM20, so that we could express the level of mitochondrial membrane potential per unit of mitochondrial membrane area measured. For a detailed description (see Supplementary Material).

Prior to assessing statistical significance, and in order to identify the appropriate statistical test, we tested the normal distribution of the sample population using the Shapiro–Wilk and Kolmogorov–Smirnov tests. Parametric tests were only used for datasets that passed both tests. Of note, all flow cytometry data were treated as non-normally distributed populations. Supplementary Table 2 demonstrates the statistical tests used for each assay to assess statistical significance. N represents the number of biological replicates. Each hPSC line is one biological replicate, and each round of differentiation from a given hPSC line represents one technical replicate. Data were analyzed and plotted using GraphPad-Prism (Prism 7.0, GraphPad Software, La Jolla, CA, United States).

Two hESC lines (429 and 360) and three hiPSC lines (established from two independent Detroit 551 clones (clones 7 and 10), and one CRL 2097 iPSC clone, CRL-8) were selected for mesoderm differentiation in a 96-well plate format (Balafkan et al., 2020). The differentiation process is shown schematically in Figure 1A. For comparative purposes, we divided the differentiation process into phases based on the expression of cell-type-specific markers (Figure 1A): pluripotent state (S1, day 0), mesendoderm cells (S2, day 2), cardiac mesoderm (S3, day 3), Isl1⁺ progenitor cells (S4, days 5–7), and cardiomyocyte and non-cardiac cells (S5, days 12–15). The specific markers for the mesoderm cell lineage were assigned based on previous studies (Sturzu and Wu, 2011; Vliet et al., 2012). Initial transcriptomic profiling of Isl1⁺ progenitor cells (dataset A) derived from different hPSCs was performed with samples collected from four hiPSCs and two of the hESC lines (429 and 360) at pluripotent state (S1) and Isl1⁺ progenitor cells (S4). Sample clustering based on the RNA-seq gene expression evidenced a marked transcriptional difference between S1 and S4, singling out the differentiation process as the main driver of transcriptional change (Supplementary Figure 1A). This was corroborated by the strong association between the first principal component of gene expression (explaining 36% of the variance) and the differentiation stage (linear regression model p = 71 × 10^–3). Cell lines were only weakly associated with subsequent principal components (Supplementary Figure 1B). To validate our gene expression analysis and further profile transcriptomes from critical stages of mesoderm differentiation, including S5, an independent RNA-seq experiment was performed using an hESC line (H1), which was solely used for transcriptomic analysis (dataset B) and not included in other sets of experiments. RNA samples from three independent differentiations of H1 were collected at the following stages: pluripotent state (S1), mesendoderm state (S2), cardiac progenitor (S4), and differentiated cells in S5. Both datasets revealed downregulation of pluripotent stem cell markers, while upregulated genes were enriched in pathways associated with cardiomyocyte differentiation including regulation of cardiac muscle, ventricular cardiac muscle tissue morphogenesis, sarcomere, and regulation of heart contraction (Figure 1B and Supplementary Figures 1C–E). The correct mesoderm differentiation route was confirmed using qPCR (Supplementary Figure 2A) and immunocytochemistry (Figures 1C,D). These data indicate that our model is a reliable model in which to assess mesoderm differentiation toward cardiomyocytes. Previously, we have shown that differentiating cells in a 96-well plate format results in a population of both cardiomyocytes and non-cardiomyocytes (Skelton et al., 2017; Balafkan et al., 2020). Thus, cells in S1 represent a pure population of PSCs (∼100% ± 0.6 SSEA4⁺), and those in S4 represent a relatively pure population of cardiac progenitors (∼92% ± 1 ISL1⁺) (Moretti et al., 2006), while S5 comprises two different cell populations, cardiomyocytes (20% ± 13 TNNT2⁺) and non-cardiomyocyte (TNNT2^–) cell populations (Figure 1E).

Thirteen polypeptides that are essential respiratory chain components are encoded by mtDNA. Unlike nuclear DNA, mtDNA is present in multiple copies, and its copy number can impact the levels of mitochondrial RNA transcripts available for generating respiratory chain subunits (Anderson et al., 1981). We assessed the mtDNA copy number at different stages of mesoderm differentiation using real-time PCR quantification relative to the nuclear gene APP (Tzoulis et al., 2013). This revealed a clear and progressive reduction of mtDNA copy number (≤85%) during differentiation of both hiPSCs and hESCs to the mesodermal lineage (Figure 2A and Supplementary Figure 3A). For comparative purposes, we quantified mtDNA in postmortem human heart using the same method and confirmed that the level of mtDNA in the mature tissue is at least 11–44-fold higher than in hPSCs and differentiated cells at S5 (Figure 2B).

Assessed by bulk RNA-seq, we found that the majority of genes involved in mtDNA homeostasis decreased during mesoderm differentiation, including mtDNA maintenance exonuclease 1 (MGME1), single-stranded DNA binding protein (SSBP), mitochondrial transcription factor A (TFAM), and mtDNA helicase (TWNK) (Figure 2C). We validated RNA-seq findings by using real-time PCR to quantify the expression of SSBP and POLG, which both have critical roles in mtDNA replication, and were able to show a progressive decline of SSBP with no significant change in POLG expression level (Supplementary Figures 3B–D). No evidence of qualitative damage such as mtDNA deletion was found at any stage of differentiation, which further confirms that reduction of mtDNA is not an artifact of damaged mtDNA replication machinery that can be caused, for example, by the continued culture of hPSCs (Supplementary Figures 3E,F).

The link between the TFAM protein and mtDNA has been investigated in depth (Larsson et al., 1994, 1998), and studies show that it binds mtDNA in molar quantities (Ekstrand et al., 2004; Kaufman et al., 2007; Kukat and Larsson, 2013). Using flow cytometry, we assessed TFAM both as a direct measure of a mitochondrial matrix protein and as an indirect measure of mtDNA level within each cell population co-stained with antibodies against stage-specific markers, to validate our findings from real-time PCR. We found a 58% reduction of TFAM in Isl1⁺ progenitor cells (S4) relative to the pluripotent stage (S1) (Figure 2D). The TFAM level decreased further in the TNNT2^– population but showed an increase in the TNNT2⁺ cells (Figure 2D).

Next, we examined changes in mitochondrial mass during mesoderm differentiation at the single-cell level. We used flow cytometry and co-staining with antibodies against stage-specific markers and Translocase of Outer Mitochondrial Membrane 20 (TOMM20), a well-established marker of mitochondrial mass. The TOMM20 level showed a significant fall (43%) from S1 to S4 and reached its lowest level at S5 (60%) (Figure 2E). As expected, the level of TOM20 was lower (∼40%) in non-cardiomyocytes than cardiomyocytes (TNNT2⁺), potentially indicating a lower level of mitochondria in non-cardiac cells. The fall in mitochondrial mass could also be followed by looking at VDAC expression. VDAC encodes a vital outer membrane protein that is routinely used as a mitochondrial mass marker; its expression started decreasing from S3 (Supplementary Figure 3D). Given the link between cell size and mitochondrial mass (Zheng et al., 2016), the finding of a lower mitochondrial mass was surprising; differentiated cells at S5 are much larger than PSCs (Figure 1D-iv). Interestingly, when we measured changes in mtDNA level relative to mitochondrial content, by plotting the level of TFAM (an indirect marker of mtDNA) against TOM20 (a direct marker for mitochondrial content), we found that this ratio varied during the course of mesoderm differentiation. The ratio of TFAM to TOM20 was significantly higher in differentiated cells in S5 relative to Isl1⁺ progenitor cells (S4) (Figure 2F), suggesting that the lowest level of mtDNA per unit mitochondrial mass occurred in Isl1⁺ progenitor cells and then rose in the more differentiated stage. Together, these findings suggest a progressive fall of mitochondrial content during mesoderm differentiation.

Despite analyzing a heterogeneous population (TNNT2⁺ and TNNT2^– cells), we were able to verify a significant increased expression of mtDNA genes in S5 relative to other stages of differentiation via bulk RNA-seq (Figure 3A). Demonstrating a significant increase in the expression of mtDNA genes in the differentiated cells in S5 suggests that differentiated cells rely more on mitochondrial respiration, regardless of the cell type they are committed to, and despite having much lower mitochondrial mass. Thus, we assessed mitochondrial oxidative activity by measuring OCR in undifferentiated stem cells (Figure 3C) and differentiated cells at S5, which contain beating cardiomyocytes (see Supplementary Movies 3, 4 and Figure 3D), using a Seahorse XF^e96 extracellular flux analyzer. Interestingly, despite a major drop in mitochondrial mass, no major change in basal respiration—the amount of oxygen consumed by mitochondria under basal conditions—was detected: differentiated cells at S5 137 pmol/min, undifferentiated hPSCs at S1 148 pmol/min (Figure 3E and Supplementary Figure 4A). This suggests that the remaining mitochondria have higher levels of oxygen consumption that compensates for the lower mitochondrial mass in differentiated cells. Maximal OCR and spare capacity (also known as reserve capacity) showed a slight increase in S5 compared with S1 cells (Supplementary Figures 4B,C), while the coupling efficiency—the proportion of mitochondrial respiratory chain activity used to drive ATP synthesis (Divakaruni and Brand, 2011)—was significantly higher (Figure 3F). The ratio of OCR to extracellular acidification rate (ECAR), an indicator of how much lactate is produced through glycolysis, also rose in S5 cells relative to S1, suggesting a shift to mitochondrial respiration (Figure 3G). The finding of a negative value for spare capacity in hPSC lines (360 and CRL-8) raised the possibility that uncoupling by CCCP had collapsed the membrane potential in these lines (see blue bars in Supplementary Figure 4A-iv) and that, with a respiratory chain already working maximally, there was no further reserve capacity to be used with the extra stress (Varum et al., 2011; Zhang et al., 2012b).

We also investigated mitochondrial membrane potential as a marker for mitochondrial activity. Membrane potential is generated by proton movement driven by the electron transport chain; and although mitochondrial membrane potential is an excellent marker for assessing mitochondrial function, it also reflects mitochondrial volume (Perry et al., 2011). We therefore employed a lipophilic cationic fluorescent probe, TMRM, to evaluate mitochondrial membrane potential at the single-cell level by flow cytometry; and we normalized our data to the MFI of TOM20 from the same differentiation batch to further adjust membrane potential to mitochondrial content within a single cell (see Supplementary Methods 1.10). We did not use cell stage-specific markers for this assessment, as this was done in live cells; and the result therefore reflects the membrane potential of total live cells in S1 and S5. The absolute median MFI of TMRM was significantly lower at S5 relative to S1 (Figure 3H and Supplementary Figure 4E); however, when we adjusted the MFI of TMRM for mitochondrial mass, using the ratio of MFI-TMRM to MFI-TOM20, we found a more than twofold increase in mitochondrial membrane potential at S5 compared with undifferentiated hPSCs at S1 (Figure 3I and Supplementary Figure 4F). The significantly higher mitochondrial membrane potential per unit of mitochondrial mass in differentiated cells at S5 relative to undifferentiated hPSCs further suggested that the mitochondria in differentiated cells are more efficient in generating ATP.

We examined cells collected on S1, S4, and S5 by TEM to assess the mitochondrial morphology and crista structure. The mitochondria in S1 appear to have wide cristae with a dense matrix, while in S4 and S5, cristae appear more compact with a clearer matrix (Figure 3J). We also identified a significant increase in the expression of genes associated with mitochondrial biogenesis and respiration activities (e.g., PPARA, PPARG, PGC-1A, and ESSRA; see Figure 3B and Supplementary Figure 5). Since no increase in mitochondrial mass was demonstrated relative to Isl1⁺ progenitor cells (Figure 2E), together, these data suggested remodeling of mitochondrial membrane and cristae as a potential explanation for the increase in mitochondrial respiration.