Comparing the signaling and transcriptome profiling landscapes of human iPSC-derived and primary rat neonatal cardiomyocytes

The inaccessibility of human cardiomyocytes significantly hindered years of cardiovascular research efforts. To overcome these limitations, non-human cell sources were used as proxies to study heart function and associated diseases. Rodent models became increasingly acceptable surrogates to model the human heart either in vivo or through in vitro cultures. More recently, due to concerns regarding animal to human translation, including cross-species differences, the use of human iPSC-derived cardiomyocytes presented a renewed opportunity. Here, we conducted a comparative study, assessing cellular signaling through cardiac G protein-coupled receptors (GPCRs) in rat neonatal cardiomyocytes (RNCMs) and human induced pluripotent stem cell-derived cardiomyocytes. Genetically encoded biosensors were used to explore GPCR-mediated nuclear protein kinase A (PKA) and extracellular signal-regulated kinase 1/ 2 (ERK1/2) activities in both cardiomyocyte populations. To increase data granularity, a single-cell analytical approach was conducted. Using automated high content microscopy, our analyses of nuclear PKA and ERK1/2 signaling revealed distinct response clusters in rat and human cardiomyocytes. In line with this, bulk RNA-seq revealed key differences in the expression patterns of GPCRs, G proteins and downstream effector expression levels. Our study demonstrates that human stem cell-derived models of the cardiomyocyte offer distinct advantages for understanding cellular signaling in the heart.

Generally, the study of human cardiac physiology and pathophysiology has been challenging. The most notable barrier has been the inaccessibility of human cardiac tissue and cell sources for in vitro experimental studies. Human heart tissue obtained from diseased hearts only provides insight into end-stage heart disease and is generally not informative of disease evolution through time. The few adult primary human cardiomyocytes that may be rarely obtained, only last a few days in culture, thereby limiting the types of experiments that can be conducted and reducing the chance to explore disease progression or therapeutic intervention. Recently, an isolation protocol has been described that was shown to maintain human primary cardiomyocytes in culture for 7 days as well as a procedure for successful cryopreservation and thawing 1 . Needless to say, the small quantities obtained still represent a substantial bottleneck beyond the reach of most laboratories.
No model is perfect. To model the human heart, other representative cell sources were exploited, most prominently being cardiomyocyte cultures derived from rodents. With large litters and relatively low-cost for maintenance compared to large-animal models, the isolation of cardiomyocytes from these species provides enough functional cells for study. Frequently, neonatal rat cardiomyocyte cultures are used as they can typically be kept in culture longer compared to adult rodent cultures that tend to de-differentiate after 2-3 days in vitro. Numerous reports have demonstrated that cardiomyocytes isolated from mice and rat hearts are signaling competent and exhibit hypertrophy in response to known inducers 2,3 . With many genetic and phenotypic similarities to humans, rodents have helped researchers study numerous aspects of cardiac physiology. However, rodent hearts do not completely recapitulate conduction and electrophysiological properties of the adult human heart. Most

Distinct protein kinase activation patterns observed in hiPSC-CM and RNCM cultures. Nuclear protein kinase A (PKA) signaling.
To compare the cellular contexts of hiPSC-CMs against rat neonatal cardiomyocytes, we began probing for differences in cellular signaling. Here, we set out to explore how similar RNCMs and immature hiPSC-CMs were. With this in mind, we measured protein kinase A activation in both cell types, at the single cell level, every 10 min for 70-min, as in the experimental pipeline depicted in Fig. 1A. To probe agonist-induced PKA activation, cardiomyocyte cultures were transduced with AAV2/6-ExRai-AKAR2-NLS (Fig. 1B,C). Three days later, cardiomyocyte cultures were imaged using an automated high content imaging system. Baseline measurements were recorded prior to stimulating cells with a panel of saturating doses of GPCR agonists targeting either endothelinergic, adrenergic or angiotensinergic systems and their downstream effectors. (Fig. 1D). As seen in Fig. 1E, the expression of the nuclear-localized PKA biosensor, as represented on the x-axis, was considerably higher in hiPSC-CMs compared to RNCMs, despite being transduced for the same length of time. When testing the relationship between biosensor expression (relative fluorescence units, RFU) and biosensor output (%ΔF/F), we did not observe a significant linear correlation in the RNCM cultures as shown with small R correlation coefficient and nonsignificant p-value (p = 0.23). In the hiPSC-CMs, the R coefficient was small, yet the p-value was significant, an observation that suggests greater variability around the regression line and a larger prediction interval when correlating both variables. Thus, to have a more uniform population of RNCMs and hiPSC-CMs, removing plausible confounds caused by biosensor expression, we set a threshold and analysed CMs on the lower end of the spectrum. Cardiomyocytes that expressed the PKA biosensor between 0 and 5000 RFU were carried forward for further analysis, and no cells exhibited zero intensity as the lowest unit measured was 462 RFU (Fig. 1E, inset).
Next, we averaged the nuclear PKA responses of all CMs, across our biological replicates, and plotted responses as a function of time (Supp. Fig. 1). RNCMs and hiPSC-CMs displayed distinct ensemble PKA signatures when stimulated by various agonists, when evaluating both kinetics and magnitude of the responses. In both CM types, forskolin, a direct activator of adenylyl cyclase, was the most potent activator of the PKA pathway followed by epinephrine. Markedly, epinephrine and norepinephrine led to larger nuclear PKA responses in RNCMs. Phenylephrine and ET-1 also resulted in PKA activation; an observation examined in detail later. To move beyond averaged responses, and as means to fully appreciate subtle differences or granularity between RNCMs and hiPSC-CMs, we used a single cell analytical approach to our biosensor datasets 10 . To identify cardiomyocyte populations or clusters that responded differently within the overall population, we merged RNCM and hiPSC-CM datasets and ran them through a clustering algorithm 11 . This revealed that the overall population of cardiomyocytes could be sub-divided into 4 main groups (Supp. Fig. 2). We arbitrarily named them, 'sustained-responders' , 'transient-responders' , 'non-responders' and 'negative-responders' . To test whether the expression of the biosensor influenced clustering, cardiomyocytes at different fluorescent intensity thresholds were independently ran through the algorithm (Supp. Fig. 2A). The inclusion of cells at varying fluorophore intensities did not appreciably impact response clusters. We also noticed that while forcing the clustering algorithm to generate four clusters instead of three (Supp. Fig. 2B), the data was partitioned into more informative patterns. By using four clusters, we noted that the population of 'responders' could be further sub-divided based on response kinetics, with some cells exhibiting 'transient' responses, while others displayed 'sustained' effects. Four clusters were carried forward in further analyses. As another internal control, we investigated whether all cardiomyocytes seeded in the microwell were exposed to a uniform concentration of drug. When stimulating the cells, we were careful to break the surface tension barrier with the microtip to ensure effective drug diffusion, as we were working with small, 10 μL volumes. The 21 fields imaged were plotted against the percentage change www.nature.com/scientificreports/ in response compared to baseline, %ΔF/F, post-drug stimulation (Suppl. Fig. 3). No significant differences were observed in response to stimulation with norepinephrine while responses to forskolin and isoproterenol showed minor positional effects by field. Still, these observations may likely be correlated with the proportion of cells that fall within a response cluster, irrespective of drug diffusion. Thus, the fraction of cells exhibiting distinct behaviors, categorized as non-responders or responders may be skewing the average response per field. To visualize the overall spread of single nuclei PKA responses, we generated heatmaps with single nuclei responses plotted on the y-axis as a function of time (x-axis). In hiPSC-CMs, adrenergic targeting ligands, isoproterenol, norepinephrine, and epinephrine, resulted in modest activation of nuclear PKA ( Fig. 2A). More specifically, the distribution revealed that 18.72 ± 5.73% of hiPSC-CMs exhibited a sustained response to isoproterenol compared to 2.66 ± 1.5% of RNCMs (Fig. 2B). Further, in response to epinephrine, 57.83 ± 2.3% of hiPSC-CMs showed a transient response compared to 35.51 ± 13.5% of RNCMs. Forskolin elicited a similar response profile in both cell types with slightly more sustained responders in the RNCM population. Most remarkably, phenylephrine resulted in nuclear PKA activity in RNCMs but not in hiPSC-CMs (Fig. 2B). A greater density of α 1 AR in RNCMs could explain the prominent PKA response upon stimulation with phenylephrine. Similarly, ET-1 drove transient PKA activation in RNCMs and hiPSC-CM, albeit to a smaller degree in hiPSC-CMs. As the ET A receptor is canonically coupled to Gα q , this observation was unexpected. However, this was previously observed in both HeLa cells overexpressing the ET A as well as in rat aortic smooth muscle cells 12 . In these models, PKA activation was reported to be independent of cAMP production. Further, in rat aortic vascular smooth muscle cells, ET-1 coupling to Gα i was shown to result in transient PKA activation, with a return to baseline after 20 min 13 . An observation analogous to what we observed in RNCM cultures with approximately 45 ± 15% of cells exhibiting this transient behaviour (Fig. 2B). G proteins may be expressed at distinct levels in both cell types, examined further below. This observation contrasts with our previous report which demonstrated phenylephrine but not ET-1 was able to induce nuclear PKA responses in RNCMs 14 . However, as ExRai-AKAR2-NLS is more sensitive than the AKAR2-NLS biosensors used in our previous study, this discrepancy may simply reflect differential sensitivities of the biosensors themselves. With a greater proportion of hiPSC-CMs exhibiting a decrease in PKA activity compared to baseline, this behaviour may suggest that adrenergic receptors exhibit a greater www.nature.com/scientificreports/ degree of constitutive activity in these cells and that in some cells, receptor occupancy decreases this. Lastly and as expected, agonists for the AT 1 R did not result in PKA activation (Suppl. Fig. 4). Predictably, Ang II and related compounds led to negligible responses in hiPSC-CMs. PMA, a direct activator of PKC, also failed to trigger PKA activation, except for a small population of hiPSC-CM. This may reflect signaling crosstalk 15 . Altogether, these trends indicate distinct signaling landscapes in rat versus human cardiomyocytes.
Nuclear ERK 1/2 signaling. We were next interested in probing for activation of the mitogen-activated protein (MAP) kinase, ERK 1/2 . Correspondingly, RNCMs and hiPSC-CMs were transduced with AAV2/6-EKAR-NLS to measure nuclear ERK 1/2 activity (Suppl. Fig. 5A,B). As observed in our PKA experiments, hiPSC-CMs expressed the EKAR-NLS biosensor at higher donor (CFP) intensities compared to RNCMs (Suppl. Fig. 5C). We filtered our CM populations using those that expressed the biosensor between 2000 and 5000 RFU. The response clusters did not appear to change greatly over this range of CFP intensities. However, before dissecting single cell response profiles, we plotted averaged ERK 1/2 responses over time. Our data demonstrated that PMA and ET-1 were the most robust drivers of nuclear ERK 1/2 in both RNCMs and hiPSC-CMs, displaying the greatest response magnitudes, while epinephrine, norepinephrine and phenylephrine exhibited significantly larger ERK 1/2 activation in RNCMs compared to hiPSC-CMs (Suppl. Fig. 5D). Next, we plotted single cell ERK 1/2 responses (Figs. 3, 4, Suppl. Fig. 6), heatmaps shown in Suppl. Figs. 7 and 8. In hiPSC-CM profiles, adrenergic ligands displayed considerable variability, ~ 30% of nuclei exhibited a decrease in ERK 1/2 activation compared to baseline. There may again be a greater degree of constitutive receptor activity in this context, perhaps via a sub-population of cells signaling through a β 2 AR-Gα i pathway. In hiPSC-CMs, ET-1 treatment led to the most robust and sustained response over time in 97.8 ± 0.5% of the nuclei population. However, even vehicle showed ERK 1/2 activity with considerably smaller magnitude suggesting there in a tone to ERK 1/2 signaling even in the absence of ligands. To further dissect this response range, we applied the clustering Single nuclear PKA data summarized as %ΔF/F (y-axis-each line is a different cell) as a function of time (x-axis). The data was partitioned into four clusters representing distinct nuclear behaviors. Nuclei exhibited either sustained or transient responses to agonists while other nuclei failed to respond or experienced a decrease in activity compared to baseline. Representation of the four response clusters was plotted as a bar chart with percentages of nuclei belonging to each response cluster as observed in (C) hiPSC-CMs and (D) RNCMs. Experiments were performed using RNCMs isolated from 4 different neonatal rat pup litters and 3 independent cardiomyocyte differentiations from hiPSCs. We observed a significant difference between species (p < 0.05) when the hiPSC-CMs or RNCMs were treated with phenylephrine. There is a significant population of cells that experienced a 'negative' response in the hiPSC-CMs (p < 0.05). www.nature.com/scientificreports/ algorithm on a select population, those previously categorized as 'responders' . This yielded two more sub-clusters, arbitrarily named 'low-responders' and 'high-responders' (Suppl. Fig. 6C). This extra step allowed us to present the data in a way that better reflected the spread as observed in the heatmaps. In hiPSC-CMs, phenylephrine effects were less pronounced and not significantly higher than vehicle (Fig. 3A, Suppl. Fig. 7). In contrast, other adrenergic ligands resulted in significant nuclear ERK 1/2 activity, specifically in response to epinephrine (p = 0.033) and norepinephrine (p = 0.038) in RNCMs compared to hiPSC-CMs (Fig. 3B, Suppl. Fig. 7). In RNCMs, transient activation patterns were observed in response to isoproterenol (p = 0.047), an observation that was not mirrored in hiPSC-CMs. Phenylephrine drove ERK 1/2 activation to a greater extent in RNCMs compared to hiPSC-CMs, once again suggesting that hiPSC-CM and RNCMs may express α 1 ARs at different levels. In adult mouse cardiomyocyte cultures, phenylephrine was shown to result in ERK activation in 60% of CMs, as α 1 AR receptors are expressed in 60% of ventricular cardiomyocytes 16,17 . This observation is similar to our ERK 1/2 assays conducted in hiPSC-CMs, at 58% (Fig. 3A). This finding is aligned with reports suggesting that mouse (rather than rat) CM cultures are better mimics of the human CM, when studying α 1 AR-mediated myocardial effects, as α 1 AR expression levels are more comparable. In contrast, rats express 5-10-fold more α 1 AR, a finding that is aligned with the strong phenylephrine responses observed in our RNCM experiments 18,19 (Fig. 3B).
Unsurprisingly, forskolin failed to produce significant ERK 1/2 signals in either CM model while phorbol 12-myristate 13-acetate (PMA) drove activation in 99% of all CMs. Similarly, Ang II and the β-arrestin biased drugs, SI and SII did not produce significant nuclear ERK 1/2 activity in comparison to vehicle (Fig. 4). It is possible that by the time our cells were imaged at the 10-min time point, the early ERK 1/2 activation phase had passed, as western blots showed ERK 1/2 activity as early as 5 min in RNCMs (Suppl. Fig. 9). However, in the abovementioned single cell FRET study, the authors also reported the lack of an Ang II-mediated ERK response 16,17 . Bulk RNA-seq reveals transcriptome-level differences in three cell types: hiPSC-CMs, RNCMs and HEK 293 cells. GPCR-mediated signal transduction-related gene sets. As distinct signaling signatures were observed above, we next sought to test whether these differences could be explained via differential gene The three response clusters were plotted as stacked bar charts with percentages of nuclei belonging to each cluster as observed in (A) hiPSC-CMs and (B) RNCMs. Insets reflect sub-clustering applied to the 'responding' nuclei cluster. It is apparent that within the 'responding' sub-population, two other clusters can be identified based on their responding magnitudes-low or high responders. Experiments were performed using RNCMs isolated from 3 different neonatal rat pup litters and 3 independent cardiomyocyte differentiations from hiPSCs. We observed a significant difference between species (p < 0.05) when the hiPSC-CMs or RNCMs were treated with isoproterenol, in both the responders and non-responder's category. The same was observed when the cells were treated with epinephrine and norepinephrine. www.nature.com/scientificreports/ expression. Thus, we set out to determine if RNCM and hiPSC-CMs were wired differently, at the endogenous level. To characterize the transcriptional landscapes in hiPSC-CM and RNCM function, we surveyed genes linked with GPCR signaling at baseline, independent of ligand. Since HEK 293 cells have been a valuable model system for dissecting GPCR function in vitro, we included these cells in our analysis 20 . Basal abundance of select cardiac-relevant Class A GPCRs, reported as transcript per kilobase million (TPM), revealed striking differences between the three in vitro cell models assessed (Fig. 5). As suggested earlier, indeed, the α 1 -adrenergic receptor was more highly expressed in RNCM compared to hiPSC-CMs or HEK 293 cells (Fig. 5A, Supp. Fig. 10A). The β 1 -adrenergic receptor was also more abundant in RNCMs. The endothelin A (ET A ) receptor was expressed at higher levels in hiPSC-CMs compared to RNCMs. In contrast, angiotensin II type I (AT 1 ) receptors were expressed at low levels in all three cell types. Certain aspects of cardiac remodeling associated with the angiotensin system have been reported to occur through fibroblast activation and subsequent paracrine signaling 21 . Regarding the endogenous signaling machinery, Gα s was shown to be well expressed in all three cell types (Fig. 5B, Supp. Fig. 10C). Gα q was expressed at comparable levels in all three cell types, though at slightly higher levels in RNCMs. Gα i2 expression threefold higher in RNCMs. Gβ isoforms were all observed, albeit to varying degrees, with Gβ1 being most abundant followed by Gβ2 (Fig. 5C, Supp. Fig. 10D). Equally, all Gγ isoforms were detected, again to different levels, with Gγ5, 10 and 12 being most abundant (Fig. 5D, Supp.  Fig. 10E). Lastly, distinct expression patterns were observed for select GPCR effector molecules (Fig. 5E, Supp.  Fig. 10F). For example, adenylyl cyclase isoforms 5 and 6 were more abundant in cardiomyocytes compared to HEK 293 cells, an observation consistent with previous reports 22 . In contrast, β-arrestin2 was expressed at greater levels in HEK 293 cells, perhaps affecting GPCR desensitization or post G protein signaling. PKA subunits were well represented, as were members of the MAPK cascade, albeit at different levels (Fig. 5E, Supp. Fig. 10F). All in all, mRNA levels seem to align well with the cellular signaling signatures reported above. Further, differential DESeq2 analysis demonstrated that genes that were either significantly up-or down regulated in RNCMs and hiPSC-CMs after 60-min agonist stimulation differed (Suppl. Fig. 11). A 60-min time point was simply selected to parallel the time course of our signaling experiments. This finding indicates that their transcriptional landscapes are also different, driving distinct cellular effects and possibly fates. Taken together, evidence presented implies www.nature.com/scientificreports/ that cell context is a critical determinant dictating receptor outcomes. For physiologically relevant conclusions to be drawn, it is important to consider the cell context and its influence especially when translating research findings for human predictability 23 .
Cardiomyocyte-, maturity-and metabolism-related gene sets in hiPSC and RNCMs. We next examined gene sets associated with cardiomyocyte fate, maturity and 'adult-like' identity. When observing normalized transcript abundance of a curated gene list for cardiac progenitor differentiation, RNCMs as well as hiPSC-CMs expressed GATA4, Nkx2-5 and Mef2c, key transcription factors that control cardiomyocyte fate 24,25 (Suppl. Fig. 12A). Genes encoding sarcomeric proteins were also abundant including TNNT2, TNNi3, Actc1 and MYH6. In our cultures, the fetal troponin isoform, TNNI1 was more abundant in hiPSC-CMs compared to the adult equivalent, TNNI3 (Suppl. Fig. 12B). RNCMs express approximately fivefold more TTNI3 and 2.5-fold more TNNI1 compared to day 28 hiPSC-CMs. This data is consistent with literature suggesting that hiPSC-CMs have a more neonatal as opposed to adult cardiomyocyte character. Both hiPSC-CMs and RNCMs express MYH6, present in the developing ventricle and neither expressed the adult isoform MYH7. This echoes the neonatal nature of both these cardiomyocyte models. Genes involved in glycolytic metabolism, such as ALDOA, BPGM, ENO1 and LDHA were more abundant in RNCM compared to hiPSC-CMs (Suppl. Fig. 13A,B) 26 . These genes have been shown to be downregulated in 'matured' hiPSC-CMs that undergo a metabolic switch from glycolysis to fatty acid oxidation 27,28 . Overall, transcript reads were higher in glycolytic metabolism pathways compared to gene sets of fatty acid related metabolism. This result is suggestive that both RNCM and hiPSC-CMs use glycolysis as an energy source.
In the context of cardiomyocyte functionality and calcium regulation, both RNCMs and hiPSC-CMs express key genes: RyR2, Cacna1c, Camk2d, Atp2a2 (SERCA2), Casq2, Slc8a1 (NCX1) and Tmp1 (Fig. 6A). Transcript levels for Bin1, a gene that participates in T-tubule formation was low in both RNCM and hiPSC-CMs, an observation consistent with the literature (Fig. 6A,B). In connection with genes involved in propagating the cardiac action potential, the voltage-gated potassium channel, KCNH2 (hERG), was more abundant in hiPSC-CMs www.nature.com/scientificreports/ than RNCMs (Fig. 6C, Suppl. Fig. 14). The HCN4 subunit of the funny current, I f , was shown to be expressed at twofold higher abundance in hiPSC-CMs compared to RNCMs. This observation is in line with the automaticity seen in hiPSC-CMs cultures and not in RNCMs 29 . The low (negligible: 0.35) transcript levels of KCNJ2 encoding for the inward rectifier potassium current I K1 has also been associated with the spontaneous contractile properties of hiPSC-CMs 30 . In summary, hiPSC-CMs share common and distinct features to RNCMs but are similar in terms of maturity.

Discussion
Translational research programs take advantage of numerous immortalized, primary cell and now hiPSC-derived cellular models to mimic conditions relevant to human disease 31 . Several reports have attested to the usefulness of hiPSCs and their differentiated derivatives as human model systems especially when attempting to recapitulate disease in a dish 32,33 . Still, we thought it important to perform a comparative study with primary rat cardiomyocytes as hiPSC-CMs are associated with certain limitations, i.e., their inherent immaturity 34,35 . As such, we combined single cell signaling with bulk transcriptomics to compare the traditional rat neonatal cardiomyocyte with human iPSC-CMs. Overall, our data serves to suggest that RNCMs and hiPSC-CMs exhibit unique signaling signatures that are, in part, rationalized by cross-species transcriptomic differences. These differences were captured in non-matured hiPSC-CMs, so that even if they are not adult-like, they are still superior human proxies compared to traditional rat primary cardiomyocytes. Selecting the correct model system becomes even more relevant in the context of cardiovascular disease as receptor and effector expression levels have been reported to be altered. In the human failing heart, for example, the ratios of β-adrenergic receptors, G proteins, and GRKs are different compared to non-failing hearts 36 . These changes affect cellular signaling profiles. Therefore, working with hiPSC-CMs enables patient access, and allows a glimpse at how cellular context affects signaling in disease. As hiPSC-CMs endogenously express the main www.nature.com/scientificreports/ therapeutic targets in heart failure, there is no need for overexpression of recombinant proteins, as commonly done in HEK 293 cells. Even if RNCMs express the relevant players and effectors, their levels of expression and stoichiometries are different. Since their limited time in culture poses an extra experimental constraint, hiPSC-CMs represent an excellent proxy for modeling the human healthy or diseased heart in an in vitro setting. hiPSCs are also amenable to repeated measures, permitting experiments that explore the effects of prolonged or repeated drug exposures on numerous cellular phenotypes. This study has limitations. It is important to consider certain confounding variables when analysing our single cell signaling datasets. The rat neonatal cultures were essentially in an in vivo context 5 days prior to the assay and data collection while hiPSC-CMs were in an in vitro context for 40+ days before being assayed. Secondly, even if the RNCM maintenance media was supplemented with a mitotic inhibitor, one must be aware that a small population of cardiac fibroblasts renders the overall population somewhat impure. However, our culture method has been shown to yield a population composed of ≥ 90% cardiomyocytes. Also, the differentiation protocol used in this study yields a mixed population of ventricular, atrial and pacemaker cardiomyocytes 8 . Thus, some of the clusters observed may reflect these different cell types exhibiting distinct behaviours [37][38][39] . A single cell sequencing experiment would allow us to further resolve the identity of the clusters in a subsequent study. Alternatively, it may be possible to use lineage-specific promoters or post-hoc analyses using antibodies. Thirdly, it would be interesting to assess how matured hiPSC-CMs compare to those used here as well as to RNCMs or adult primary cultures 40 . Nevertheless, the main message of our study remains, that the correct stoichiometry of receptors and effectors is critical for meaningful conclusions to be drawn from the data particularly for translational efforts. All in all, our study demonstrates that human stem cell-derived models of the cardiomyocyte do provide significant advantages, even if immature, and should be taken advantage of if their use is appropriate and aligned with research objectives.

Materials and Methods
Reagents. Unless specified, all common laboratory reagents were purchased from Sigma-Aldrich. Drugs Rat neonatal cardiomyocyte isolation and culture. All procedures involving animals were approved by the McGill University Animal Care Committee, in accordance with Canadian Council on Animal Care Guidelines. Sprague-Dawley dams with postnatal day 1-3 pups were purchased from Charles River, Saint-Constant QC, Canada. Complete litters, including both female and male neonatal pups were sacrificed by decapitation, as previously described 14,41,42 . Using forceps, whole hearts were isolated from the chest cavity and placed in cold Hank's balanced salat solution (HBSS, Wisent, 311-511-CL). Using surgical scissors, hearts were kept whole and cut 3-5 times to increase the surface area for overnight enzymatic digestion at 4 °C with 0.1% trypsin in HBSS. The next morning, the trypsin reaction was inhibited with 7% FBS supplemented DMEM low glucose + penicillin/streptomycin (P/S). Five serial collagenase digestions were then performed. Cell suspensions containing cardiomyocytes and cardiac fibroblasts from whole heart extracts were then seeded onto tissue-culture treated plastic 10-cm dishes. Non-cardiomyocyte cells such as fibroblasts can attach to plastic while cardiomyocytes do not adhere to plastic surfaces. The cell suspension after this incubation period is thus enriched for cardiomyocytes. After two 75-min incubations, the cell suspension typically contains > 90% cardiomyocytes and were seeded in black optical bottom 96-well plates (Thermo Scientific, 165305) coated with human plasma fibronectin + 0.1% gelatin in DMEM low glucose + 7% (vol/vol) FBS + P/S + 10 μM cytosine-β-d-arabinoside (AraC, Sigma Aldrich, C1768-500MG). The mitotic inhibitor, AraC, was added to prevent proliferation of remaining fibroblasts. Cardiomyocytes were then maintained in a humidified atmosphere at 37 °C with 5% CO 2 . The next day, cardiomyocytes were washed three times with DMEM low glucose and exchanged with maintenance medium containing DMEM low glucose and universal ITS supplement composed of insulin, transferrin, and selenium (Wisent, 315-081-QL) prior to being transduced with adeno-associated viruses carrying the genetically encoded biosensor of interest. For the RNAseq experiment, RNCMs were seeded in a 6-well dish coated with human plasma fibronectin + 0.1% gelatin. After vehicle or agonist stimulation, RNCMs were incubated for 1 h at physiological conditions whereafter RNCMs were harvested, and RNA was isolated.

Differentiation of human induced pluripotent stem cells into cardiomyocytes. The use of human derived iPSCs in this research was approved by the McGill University Health Centre Research Ethics
Board. Here, an hiPSC line derived from a healthy donor (AIW002-2) was provided from the Montreal Neurological Institute through the Open Biorepository, C-BIGR 43,44 . Human iPSCs, between passages 2 to 8 post-thaw, were differentiated into cardiomyocytes following the established GiWi protocol with slight modifications 45 . hiPSC-CMs were routinely screened for mycoplasma contamination and both hiPSCs and differentiated hiPSC-CMs used in this study were mycoplasma-free. For cardiomyocyte lineage commitment, a monolayer of 500,000 hiPSCs, dissociated with Accutase, were seeded onto Matrigel (Corning, 354277) coated 24-well dishes in mTeSRPlus media supplemented with Rho kinase inhibitor, Y-27632-HCl (Selleckchem, S1049). The next day, the media was exchanged for fresh mTeSRPlus media. On day 0 of the protocol, iPSCs were stimulated with Wnt activator, CHIR99021 (Cayman Chemical, 13122) for 24 h in RPMI 1640 media supplemented with B27 minus insulin (ThermoFisher, A1895602). At the same time the next day, the media was exchanged for fresh RPMI 1640 supplemented with B27 minus insulin. On day 3, Wnt was inhibited using IWP2 (Selleckchem, S7085). Cloning biosensors in AAV plasmids. The FRET-based ERK 1/2 protein kinase biosensor EKAR-EV was generously provided by Dr. Michiyuki Matsuda and carried a nuclear localization (NLS) signal sequence 47 . This FRET-based biosensor was introduced into an AAV compatible backbone, pENN-AAV-CAG-tdTomato (Addgene catalog #105554) using a 'cut-and-paste' method using BamHI and BstBI restriction enzymes and ligase. The ExRai-AKAR2-NLS biosensor was obtained from Dr. Jin Zhang's lab 48 . For cloning ExRai-AKAR2 biosensor into an AAV backbone, a multiple cloning site (MCS) was first generated in pAAV-CAG-hChR2-H134R-tdTomato (Addgene catalog #28017). The MCS was generated by annealing the following two primers 5′-gatccgctagcgtttaaacttaaggtaccgagctcactagtgaattctgcagatatccagcacagtggcggccgctcgagggcccttcga-3′ and 3′-gcgatcgcaaatttgaattccatggctcgagtgatcacttaagacgtctataggtcgtgtcaccgccggcgagctccc gggaagcttcga-5′. To introduce the ExRai-AKAR2 biosensor into the pAAV-CAG-MCS backbone, biosensors were amplified by PCR and inserted using 5′-NheI and 3′-HindIII. Universal forward primer sequence was 5′-gctagctagcgccaccatgctgcgtcgcgccaccctg-3′. Reverse primer was used 5′-catagaagcttttatgcgtcttccacctttc-3′ for ExRai-AKAR2-NLS. High-content imaging of primary neonatal rat and hiPSC-derived cardiomyocytes. Prior to conducting the signaling experiments, cardiomyocyte cultures were inspected under a phase-contrast microscope to ensure RNCMs and hiPSC-CMs were healthy. Visual inspection of the cardiomyocytes permitted us to confirm that the hiPSC-CMs were spontaneously contracting and that both hiPSC-CMs and RNCMs appeared to have cardiomyocyte-like morphologies. In this study, the assay buffer was clear HBSS with calcium, magnesium, and sodium bicarbonate (Wisent, 311-513-CL). Sterile assay buffer was warmed to physiological temperature, 37 °C, and cardiomyocytes were washed 2-3 times prior to imaging. Cells were left bathing in 90 μL of HBSS for imaging. Post-washing, cardiomyocyte cultures were re-incubated in a humidified atmosphere of 37 °C with 5% CO 2 for appropriately 1 h before imaging to allow cells to re-equilibrate in the assay buffer. During this time, the temperature control settings (TCO) of a Perkin Elmer Opera PHENIX high-content screening system were set allowing sufficient time for to warm up to 37 °C with 3% CO 2 for live-cell imaging. All drugs were prepared tenfold more concentrated stock as 10 μL would be added in each well representative of a tenfold dilution, 100 μL final volume. For nuclear ERK 1/2 assays, images were acquired using a 20 × air objective using a 425 nm laser for excitation of CFP. Emissions were detected with filters at 435-515 nm (CFP-donor) and 500-550 nm (YFP-acceptor). For the ExRai-AKAR2-NLS biosensor, images were acquired using a 20 × air objective using a 375 nm and 480 nm laser for excitation and 500-550 nm emission filter. Using Perkin Elmer's Harmony software, the experiment was set up in such a way where the microscope would first acquire a ligandindependent measurement before the microwell plate would be ejected allowing time to perform drug stimulations. The plate itself was not displaced during drug treatment. Post-treatment measurements were automated and continued for another 70 min for a total of 8 readings, acquired at 10-min intervals.

Transduction of primary neonatal rat and hiPSC-derived cardiomyocytes.
Image analysis and data processing in R. After live-cell imaging, images were imported into Perkin Elmer's Columbus software. In Columbus, images were processed by first identifying the nuclei of individual cardiomyocytes followed by the calculation of morphology features. This output provided information regarding individual nuclei size and shape including roundness. To remove cell debris from the analysis, a size threshold was set as well as a roundness criterion. Intensity properties were calculated for both CFP and YFP fluorophores for EKAR-NLS and the same was done for the single fluorophore ExRai-AKAR2-NLS biosensor. Next, FRET ratios were computed by dividing acceptor/donor (YFP/CFP) for EKAR-NLS and 488 nm/375 nm for the ExRai-AKAR2-NLS biosensors, respectively. Columbus was then queried to output these data for each cardiomyocyte that fit the criteria listed and data was exported as text files for further analysis. Data were then processed in R based on a previously published in-house single-cell analytical approach with slight modifications, formerly applied to neuronal cultures 10  www.nature.com/scientificreports/ ward. As nuclei tend to drift slightly overtime, occasionally due to minor microplate positional changes occurring during ligand addition, we set a parameter where the fluorescence intensity of each nucleus could only deviate by ≤ 20% intensity to be considered as the same object/nuclei. To calculate the delta FRET or change in FRET in response to drug stimulations, we subtracted the ligand induced FRET from basal, ligand independent FRET. This expression was then converted into a percentage change in FRET (%ΔF/F) where the denominator (F) represented basal FRET. For consistency, (F) was computed by averaging the basal FRET across all nuclei in the same microwell. To detect whether cardiomyocyte sub-populations exhibited differential agonist-induced responses, when comparing RNCMs and hiPSC-CMs, we applied a clustering algorithm using the TSclust package 11 in R. Nuclei responses were clustered based on their magnitude of response over time, using the 'pam' function. For the single cell clustering, hiPSC-CMs and RNCMs nuclei were merged into a single dataset. Merging both sets of data allowed for output clusters to be reliable across RNCMs and hiPSC-CMs as means to obtain matched patterns and clusters within the two populations being compared. The algorithm could yield distinct clusters when applied to independent data sets. Following clustering, the merged hiPSC-CM and RNCM datasets were split and further plotted as heatmaps for visualization using pheatmap package. Stacked bar charts were also generated to better summarize the data using the ggplot2 package.
RNAseq preparation and analysis. RNCM and hiPSC-CM cultures were seeded on fibronectin-coated 6-well dishes at a density of 1 million cells per well. After 1 h vehicle and drug stimulation, RNA was collected and isolated with the QiAshredder kit (Qiagen, 79656) followed by RNeasy Mini Kit (Qiagen, 74106) according to manufacturer's instructions. RNA quality was later assessed by using a NanoDrop to ensure for RNA purity. Libraries were prepared using the NEBNext rRNA-depleted (HMR) stranded library kit and paired-end 100 bp sequencing to a depth of 25 million reads per sample. Sequencing was performed on the Illumina NovaSeq™ 6000 at the McGill University and Genome Québec Innovation Centre. RNAseq analysis was performed as we have previously described 14,51 . Briefly, quality of reads was determined using FastQC, and trimmed with TrimGalore. Trimmed reads were then aligned to the Ensembl rat reference genome (Rattus_norvegicus.Rnor_6.0.100) or the reference human genome (Homo_sapiens.GRCh38.104) with STAR 14,51 . Transcripts were then assembled with StringTie. Normalized transcript abundance, measured as TPM, transcript per million, was plotted using pheatmap package in R. Curated gene sets were downloaded from https:// www. gsea-msigdb. org/ gsea/ index. jsp 24,25 .
Statistical analysis. All statistical analyses were performed using GraphPad Prism, where summary data processed and computed in R was imported. For cellular signaling datasets, Figs. 2, 3 and 4 and Suppl. Fig. 4, non-parametric Welch's t-tests were conducted as test statistics to determine whether RNCM and hiPSC-CM populations and clusters had equal means. In Suppl. Fig. 3, we performed Welch's ANOVA followed by Dunnett's multiple comparisons test, to test whether response means were equal across different fields measured. Welch's ANOVA was chosen as Bartlett's test statistics indicated that there was a significant difference between variances.

Animal experiments. All procedures involving animals were approved by the McGill University Animal
Care Committee, in accordance with Canadian Council on Animal Care Guidelines and in accordance with ARRIVE guidelines.
iPSCs. We have institutional approval to obtain blood and generate iPSCs and their associated differentiated cells. All enrolled participants provided informed consent with REB approval HID-B/2020-6362.

Data availability
Data available from the Gene Expression Omnibus: record GSE225707.