Single-Cell RNA-Sequencing and Optical Electrophysiology of Human Induced Pluripotent Stem Cell-Derived Cardiomyocytes Reveal Discordance Between Cardiac Subtype-Associated Gene Expression Patterns and Electrophysiological Phenotypes

Abstract

The ability to accurately phenotype cells differentiated from human induced pluripotent stem cells (hiPSCs) is essential for their application in modeling developmental and disease processes, yet also poses a particular challenge without the context of anatomical location. Our specific objective was to determine if single-cell gene expression was sufficient to predict the electrophysiology of iPSC-derived cardiac lineages, to evaluate the concordance between molecular and functional surrogate markers. To this end, we used the genetically encoded voltage indicator ArcLight to profile hundreds of hiPSC-derived cardiomyocytes (hiPSC-CMs), thus identifying patterns of electrophysiological maturation and increased prevalence of cells with atrial-like action potentials (APs) between days 11 and 42 of differentiation. To profile expression patterns of cardiomyocyte subtype-associated genes, single-cell RNA-seq was performed at days 12 and 40 after the populations were fully characterized with the high-throughput ArcLight platform. Although we could detect global gene expression changes supporting progressive differentiation, individual cellular expression patterns alone were not able to delineate the individual cardiomyocytes into atrial, ventricular, or nodal subtypes as functionally documented by electrophysiology measurements. Furthermore, our efforts to understand the distinct electrophysiological properties associated with day 12 versus day 40 hiPSC-CMs revealed that ion channel regulators SLMAP, FGF12, and FHL1 were the most significantly increased genes at day 40, categorized by electrophysiology-related gene functions. Notably, FHL1 knockdown during differentiation was sufficient to significantly modulate APs toward ventricular-like electrophysiology. Thus, our results establish the inability of subtype-associated gene expression patterns to specifically categorize hiPSC-derived cells according to their functional electrophysiology, and yet, altered FHL1 expression is able to redirect electrophysiological maturation of these developing cells. Therefore, noncanonical gene expression patterns of cardiac maturation may be sufficient to direct functional maturation of cardiomyocytes, with canonical gene expression patterns being insufficient to temporally define cardiac subtypes of in vitro differentiation.

Keywords: cardiomyocyte; differentiation; electrophysiology; gene expression; induced pluripotent stem cells.

FIG. 1.

ArcLight demonstrates the ability to detect AP morphologies consistent with current density of I_Kur. (A) ArcLight-expressing cells demonstrate decreased fluorescence with depolarization. Scale bar represents 20 μm. (B) Optical APs include those that resemble ventricular-like (D41), atrial-like (D41), and nodal-like (D42) morphologies. (C) Filtered optical tracings of negative change in fluorescence over fluorescence (−ΔF/F) allow analysis of AP properties, including APD₅₀ (horizontal green dotted line), APD₉₀ (horizontal red dotted line), amplitude (blue line), and V_max (not shown). (D) Representative traces for AP recordings before and after treatment with 200 nM DPO-1 using current clamp mode at a constant rate of 1 Hz through 5 ms depolarizing current injections of 150–300 pA. (E) Representative whole-cell outward currents before and after 200 nM DPO-1 elicited by depolarization of 300 ms duration to +40 mV from a holding potential of −50 mV. (F) Representative inverted fluorescence traces for optical APs before (black) and after (red) 200 nM DPO-1. (G) Change in APD₉₀/APD₅₀ (baseline − treatment) following DPO-1 administration, versus baseline APD₉₀/APD₅₀. Cells n = 190 from three to four differentiations each of three clones. Spearman correlation coefficient is shown. (H) Cells with baseline APD₉₀/APD₅₀ >1.4 adopt more ventricular-like AP morphology (larger change in APD₉₀/APD₅₀) following treatment with DPO-1. (I) Cells with baseline APD₉₀/APD₅₀ >1.4 exhibit shorter APD₅₀. Cells for (D–I) are between D30 and D40 of differentiation. Trace bars represent 500 ms. All data are reported as median ± interquartile range. P values calculated via Mann–Whitney U test. AP, action potential; APD₅₀, action potential duration at 50% repolarization; APD₉₀, action potential duration at 90% repolarization; D, day; I_Kur, ultrarapid delayed rectifier potassium current; V_max, maximum upstroke velocity.

FIG. 2.

hiPSC-CMs demonstrate electrophysiological maturation and increased heterogeneity with extended time in culture. (A) Optical AP amplitudes are shifted toward larger values after day 16 of differentiation. (B) Optical AP V_max is shifted toward larger values after D16 of differentiation. (C) hiPSC-CMs between D37 and D40 of differentiation (late) respond to TTX with an increased interval between APs compared with the corresponding vehicle treatment. Cells between D16 and D20 of differentiation (early) also respond, but to a lesser degree. Cells were analyzed in six independent experiments from three different clones. Data are reported as median ± interquartile range. P values between vehicle and TTX treatments were calculated via a Mann–Whitney U test. (D) AP morphology, as described by APD₉₀/APD₅₀ (ratio of action potential durations at 90% and 50% repolarization), shows a shift in distribution with increased time in culture. The horizontal line marks APD₉₀/APD₅₀ ratio of 1.4. Data for (A, B, D) were collected from four to seven independent differentiations per time range. Three clones from unrelated individuals are represented. White dots within each violin indicate medians and black rectangles indicate interquartile range. Days 11–13 of differentiation: n = 121 cells; D11–D14: n = 123; D17–D20: n = 88; D23–D27: n = 121; D32–D34: n = 102; D38–D42: n = 166. Following statistical analysis via a Kruskal–Wallis test, pairwise comparisons were performed using Dunn's test with Bonferroni adjustment. Significance of pairwise comparisons is presented as box color in the corresponding matrix for each parameter. hiPSC-CMs, human-induced pluripotent stem cell-derived cardiomyocytes; TTX, tetrodotoxin.

FIG. 3.

Single-cell RNA-seq analysis in hiPSC-derived cardiac differentiation clearly separates days 12 and 40 cells with distinct electrophysiological properties. (A) Schematic of experimental design. hiPSC-CMs were analyzed at days 12 and 40 of differentiation via single-cell RNA-seq and ArcLight. (B) Unsupervised hierarchical clustering on the entire transcriptomic profiles of all cells (D12: n = 42 cells; D40: n = 43 cells). (C) Principal component analysis compiled from the entire transcriptomic profiles of all cells. (D) Representative inverted fluorescent traces for ventricular-like morphology (from D12) and atrial-like morphology (from D40). (E) Percentage of total cells at D12 or D40 with each AP morphology classification. “Indeterminate” cells had either insufficiently mature APs or demonstrated a mixture of atrial- and ventricular-like APs. (F) AP amplitude (left), V_max (middle), and APD₉₀/APD₅₀ (right) evaluated from D12 and D40 populations before submission for sequencing. P values were calculated by either a Student's t-test or Mann–Whitney U test. Data are reported as median ± interquartile range. hiPSC, human induced pluripotent stem cell.

FIG. 4.

Single-cell RNA-seq profiles from hiPSC-CMs do not demonstrate distinct cardiomyocyte subtype commitment. (A) Heatmap of expression patterns for cardiomyocyte genes across all 85 cells on days 12 and 40 (each column representing one cell). Blue and yellow represent low and high expression, respectively. MYL2 transcripts were detected in >25% of D12 and 100% of D40 cells, NKX2.5 transcripts were detected in >75% of all cells, and all other transcripts were detected in every cell. (B) Boxed D12 and D40 cells were selected as cardiomyocytes for further analysis, according to unsupervised hierarchical clustering, robust expression of cardiomyocyte genes, and low expression of cardiac fibroblast genes. (C) Heatmap of expression patterns for atrial-, nodal-, and ventricular-associated genes in cardiomyocytes. Genes with significantly different expression (log₂FC >2 and P < 0.05) between D12 and D40 are marked with an asterisk. (D) Heatmap representation of a correlation matrix for expression of atrial- and ventricular-associated genes, labeled in green (atrial) or red (ventricular). Positive and negative correlations are labeled as orange and blue, respectively.

FIG. 5.

Several ion channel regulators exhibit increased expression with electrophysiological maturation of hiPSC-CMs. (A) Heatmap representation of expression patterns for representative electrophysiology- and calcium handling-associated genes in cardiomyocytes. Genes with significantly different expressions (log₂FC >2 and P < 0.05) between days 12 and 40 are marked with an asterisk. (B–D) Relative expression (log₂FC) of SLMAP (B), FGF12 (C), or FHL1 (D) D40 gene expression for both single-cell RNA-seq (cardiomyocytes only) and qRT-PCR, compared with D12 samples. Error bars represents SEM. qRT-PCR data were obtained from one to three clones each derived from three unrelated individuals. Student's t-test analysis was used to calculate P values associated with the qRT-PCR data. FC, fold change; qRT-PCR, quantitative reverse transcription-polymerase chain reaction; SEM, standard error of the mean.

FIG. 6.

Reduced expression of FHL1 shifts hiPSC-CM AP morphologies away from atrial-like properties. (A) Immunofluorescence of cTnT and FHL1 at day 40 (representative of three analyzed clones). Scale bar represents 50 μM. (B) Representative western blot and quantification from five clones (representing three unrelated individuals) showing FHL1 knockdown efficiency via shRNA. Data are reported as mean ± SEM. (C) Cells from clones that had a mean APD₉₀/APD₅₀ > 1.4 in the control condition had increased APD₅₀ with FHL1 knockdown. Data were collected from same differentiations as represented in (B). (D) Cells from clones that had a mean APD₉₀/APD₅₀ > 1.4 in the control condition had decreased APD₉₀/APD₅₀ with knockdown of FHL1. Data for (C, D) are reported as median ± the interquartile range. (E) Heatmap of differentially expressed genes from FHL1 knockdown with potential involvement in cardiomyocyte development (FOXP1-AS1, PROX1-AS1, WNT2, WNT2B) or function (electrophysiology: CACNA1B, CACNA1C-AS3, CACNA1C-AS4, CACNA1C-IT3; gap junction: GJA9; contraction: TTN, MYH4, TNNC2, ACTN3, MYH8, TNNT3). P values were calculated by either a Student's t-test or Mann–Whitney U test. shRNA, short hairpin RNA.