Hypertrophic cardiomyopathy-associated mutations drive stromal activation via EGFR-mediated paracrine signaling

Hypertrophic cardiomyopathy (HCM) is characterized by thickening of the left ventricular wall, diastolic dysfunction, and fibrosis, and is associated with mutations in genes encoding sarcomere proteins. While in vitro studies have used human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) to study HCM, these models have not examined the multicellular interactions involved in fibrosis. Using engineered cardiac microtissues (CMTs) composed of HCM-causing MYH7-variant hiPSC-CMs and wild-type fibroblasts, we observed cell-cell cross-talk leading to increased collagen deposition, tissue stiffening, and decreased contractility dependent on fibroblast proliferation. hiPSC-CM conditioned media and single-nucleus RNA sequencing data suggested that fibroblast proliferation is mediated by paracrine signals from MYH7-variant cardiomyocytes. Furthermore, inhibiting epidermal growth factor receptor tyrosine kinase with erlotinib hydrochloride attenuated stromal activation. Last, HCM-causing MYBPC3-variant CMTs also demonstrated increased stromal activation and reduced contractility, but with distinct characteristics. Together, these findings establish a paracrine-mediated cross-talk potentially responsible for fibrotic changes observed in HCM.

Methods
snRNA-sequencing data 

snRNA-seq dataset includes data on CMTs made with healthy wild type (WT) or hypertrophic cardiomyopathy (HCM)-causing (R403Q+/- mutation in myosin heavy chain) human induced pluripotent stem cell derived-cardiomyocytes (hiPSC-CMs) (derived from PGP1 hiPSC line) and ventricular cardiac fibroblasts (Lonza Cat. CC-2904) (1).  

Methods: A total of 60,000 cells per tissue (ten tissues pooled per sample), consisting of 90% hiPSC-CMs and 10% vCFs, were mixed in an ECM solution consisting of 4 mg/ml of human fibrinogen (Sigma), 10% Matrigel (Corning), 0.4 unit of thrombin (Sigma) per mg of fibrinogen, 5 µM Y-27632 (Tocris), and 0.033 mg/mL aprotinin (Sigma). The cell-ECM mixture was pipetted into each tissue well, and after gel polymerization for ten minutes, tissue maintenance growth media containing high glucose DMEM (Fisher) supplemented with 10% FBS (Sigma), 1% penicillin-streptomycin (P/S) (Fisher), 1% Non-essential Amino Acids (Fisher), 1% Glutamax (Fisher), 5 µM Y-27632, 0.033 mg/mL aprotinin, and 150 µg/mL L-ascorbic acid 2-phosphate sesquimagnesium salt hydrate. Y-27632 was removed two days following seeding, and the growth media was replaced every other day. 

CMTs were flash frozen in liquid nitrogen on day 7 (10 tissues pooled per sample) and stored in -80°C. Individual nuclei were isolated from frozen tissue samples. Briefly, nuclei were isolated (2) and RNA was reverse-transcribed and converted into cDNA libraries using a 10x Chromium Controller and Chromium Single Cell 3' v3.1 reagent kit (10x Genomics). Bar-coded libraries were pooled and sequenced (Illumina NovaSeq 6000). Single nucleus RNAseq alignment and gene counts were performed using Cell Ranger 1.2 (10x Genomics) and Seurat (3) and R 4.1.0 and managed via RStudio. There were 15,129 total sequenced nuclei. 

Principal component (PC) analysis was performed to determine the dimensionality of this dataset within Seurat and the number of PCs was selected using both permutation-based and heuristic methods. Cluster assignment was performed in an unsupervised manner using a shared nearest neighbor approach within Seurat. To validate the ideal number of PCs and clustering resolution in this dataset, repeat analysis was performed while varying both the dimensionality and resolution. Using the Clustree R package (4), the hierarchical evolution of cluster identity was determined. This analysis allows the user to investigate how clustering resolution impacts cell identity, in an effort to assign robust classification to cell clusters which remain consistent across resolution values. For each clustering resolution, marker genes for each cluster, TTN for cardiomyocytes and FN1 for fibroblasts, were identified in order to assign relevant biological function to various subclusters. Dimensionality and resolution were adjusted until each identified cluster contained a non-zero amount of significantly expressed marker genes that were able to assign functional gene ontological terms; in this case, at cluster resolution 0.4, which was used for clustering analysis. At greater resolution values, cluster subsets returned either too few genes to assign functional classifications or non-coding genes of unknown function. Dotplot projections of single cells were generated with UMAP coordinates, using the number of dimensions determined from PCA as described above. For the purposes of this analysis, 30 dimensions were considered for both the initial PCA, and subsequent UMAP visualization. To assign cell clusters specific identities, canonical marker gene expression was analyzed to assign broad cell classes. To assign cell subclusters, such as CM1-CM6, unbiased genetic markers for each population were calculated using the Wilcoxon Rank Sum test within Seurat, at which point at most the top 100 genes were used as a Gene Ontology query to assign cellular functions to each cell subcluster.

Genes that were upregulated R403Q+/- cardiomyocytes compared to WT cardiomyocytes with a p-value < 0.01 and fold-change greater than 1.3 were input into Enrichr (5,6).  The top upregulated pathways associated with the top upregulated genes were obtained using WikiPathways (7).

Bulk RNA-sequencing data

Bulk RNA-sequencing data include data on serum starved vCFs, vCFs treated with 100 ng/mL recombinant human epidermal growth factor (rhEGF), and conditioned media from hiPSC-CMs (1).

Methods: vCFs were trypsinized, centrifuged, and flash-frozen after two days in serum-starvation (low glucose DMEM with 0.1% fetal bovine serum) vs 100 ng/mL rhEGF or conditioned media from wild-type vs R403Q+/- hypertrophic cardiomyopathy hiPSC-CMs for three different batches. Cells were homogenized in Trizol Reagent (Life Technologies, Inc., Grand Island, NY) with TissueLyzer II (QIAGEN, Inc., Valencia, CA) and RNA was isolated by conventional methods. RNA went through two rounds of mRNA purification (polyA-selection) using Dynabeads mRNA DIRECT Kit (Invitrogen, Carlsbad, CA). Double-stranded cDNA was generated using the Superscript III First-Strand Synthesis System (Invitrogen, Carlsbad, CA). The cDNA products were used to construct libraries with the Nextera XT DNA Sample Preparation Kit (Illumina, Inc., San Diego, CA). Libraries were paired-end sequenced with a read length of 75 base pairs (75PE) on the Illumina NextSeq 500. Reads were aligned to the hg38 Human Genome using Spliced Transcripts Alignment to a Reference (STAR) (8). Data were analyzed as previously described (9) and normalized to the total number of reads per kilobase of exon per million (RPKM). 

Genes that were upregulated in rhEGF-treated vCFs compared to serum-starved vCFs with a p-value < 0.01 and fold-change greater than 1.3 were input into Enrichr (5,6). For conditioned media experiments, gene expression values were averaged over three experimental repeats. Genes that had an average fold-change expression greater than 1.15 with a p-value < 0.01 for all three experiments were input into Enrichr. The top upregulated pathways and transcription factors associated with the top upregulated genes were obtained using WikiPathways (7) and ChEA (10) databases, respectively.

1. J.K. Ewoldt et al. Hypertrophic cardiomyopathy–associated mutations drive stromal activation via EGFR-mediated paracrine signaling.Sci. Adv.10 (2024).

2. E. R. Nadelmann, J. M. Gorham, D. Reichart, D. M. Delaughter, H. Wakimoto, E. L. Lindberg, M. Litviňukova, H. Maatz, J. J. Curran, D. Ischiu Gutierrez, N. Hübner, C. E. Seidman, J. G. Seidman, Isolation of Nuclei from Mammalian Cells and Tissues for Single-Nucleus Molecular Profiling. Curr Protoc 1(2021).

3. Y. Hao, S. Hao, E. Andersen-Nissen, W. M. Mauck, S. Zheng, A. Butler, M. J. Lee, A. J. Wilk, C. Darby, M. Zager, P. Hoffman, M. Stoeckius, E. Papalexi, E. P. Mimitou, J. Jain, A. Srivastava, T. Stuart, L. M. Fleming, B. Yeung, A. J. Rogers, J. M. McElrath, C. A. Blish, R. Gottardo, P. Smibert, R. Satija, Integrated analysis of multimodal single-cell data. Cell 184, 3573-3587.e29 (2021).

4. L. Zappia, A. Oshlack, Clustering trees: a visualization for evaluating clusterings at multiple resolutions. Gigascience 7 (2018).

5. E. Y. Chen, C. M. Tan, Y. Kou, Q. Duan, Z. Wang, G. V. Meirelles, N. R. Clark, A. Ma’ayan, Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics 14, 128 (2013).

6. Z. Xie, A. Bailey, M. V. Kuleshov, D. J. B. Clarke, J. E. Evangelista, S. L. Jenkins, A. Lachmann, M. L. Wojciechowicz, E. Kropiwnicki, K. M. Jagodnik, M. Jeon, A. Ma’ayan, Gene set knowledge discovery with Enrichr. Curr. Protoc. 1, e90 (2021).

7. A. R. Pico, T. Kelder, M. P. Van Iersel, K. Hanspers, B. R. Conklin, C. Evelo, WikiPathways: Pathway editing for the people. PLOS Biol6, 1403–1407 (2008).

8. A. Dobin, C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski, S. Jha, P. Batut, M. Chaisson, T. R. Gingeras, STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

9. D. C. Christodoulou, H. Wakimoto, K. Onoue, S. Eminaga, J. M. Gorham, S. R. DePalma, D. S. Herman, P. Teekakirikul, D. A. Conner, D. M. McKean, A. A. Domenighetti, A. Aboukhalil, S. Chang, G. Srivastava, B. McDonough, P. L. de Jager, J. Chen, M. L. Bulyk, J. D. Muehlschlegel, C. E. Seidman, J. G. Seidman, 5’RNA-Seq identifies Fhl1 as a genetic modifier in cardiomyopathy. Journal of Clinical Investigation 124, 1364–1370 (2014).

10. A. Lachmann, H. Xu, J. Krishnan, S. I. Berger, A. R. Mazloom, A. Ma’ayan, ChEA: Transcription factor regulation inferred from integrating genome-wide ChIP-X experiments. Bioinformatics 26, 2438–2444 (2010).