Single Nucleus RNA-sequencing Reveals Altered Intercellular Communication and Dendritic Cell Activation in Nonobstructive Hypertrophic Cardiomyopathy

End stage, nonobstructive hypertrophic cardiomyopathy (HCM) is an intractable condition with no disease-specific therapies. To gain insights into the pathogenesis of nonobstructive HCM, we performed single nucleus RNA-sequencing (snRNA-seq) on human HCM hearts explanted at the time of cardiac transplantation and organ donor hearts serving as controls. Differential gene expression analysis revealed 64 differentially expressed genes linked to specific cell types and molecular functions. Analysis of ligand-receptor pair gene expression to delineate potential intercellular communication revealed significant reductions in expressed ligand-receptor pairs likely affecting the extracellular matrix, growth factor binding, peptidase regulator activity, platelet-derived growth factor binding and protease binding in the HCM tissue. Changes in Integrin-β1 receptor expression were responsible for many observed changes related to extracellular matrix interactions, by increasing in dendritic, smooth muscle and pericyte cells while decreasing in endothelial and fibroblast cells, suggesting potential mechanisms for fibrosis and microvascular disease in HCM and a potential role for dendritic cells. In contrast, there was an increase in ligand-receptor pair expression associated with adenylate cyclase binding, calcium channel molecular functions, channel inhibitor activity, ion channel inhibitor activity, phosphatase activator activity, protein kinase activator activity and titin binding, suggesting important shifts in various signaling cascades in nonobstructive, end stage HCM.


Introduction
Hypertrophic Cardiomyopathy (HCM) is an inherited disorder characterized by unexplained left ventricular hypertrophy, often asymmetric, involving the interventricular septum, associated with left ventricular outflow tract (LVOT) obstruction, fibrosis, microvascular occlusion, and sudden cardiac death. Genetic studies have identified numerous causal mutations in a variety of sarcomere genes such as MYH7, MYL2, MYL3, MYBPC3, TNNT2, TNNI3 and TPM1, leading to the concept that HCM is a disease of the sarcomere [1], but genetic testing is only informative in approximately 30% of probands, indicating that other genes and factors significantly contribute to the HCM phenotype. Recently, it has been recognized that in most cases HCM can be considered polygenic with multiple loci contributing to the phenotype or in some cases acting as modifiers of existing sarcomere mutations, by affecting modifiable risk factors such as diastolic blood pressure [2][3][4]. Patients with sarcomere gene mutations have been found to have more adverse events than those without sarcomere mutations [5], thus implicating a need for better understanding of distinct pathogenic steps in these two groups. The development of fibrosis is associated with an increased risk of sudden cardiac death and is a poor prognostic factor [6][7][8][9][10][11]. Excess deposition of extracellular matrix (ECM) proteins can cause cardiac stiffening, which can impede contraction [6]. It can additionally disrupt electrical coupling, heightening patient susceptibility to arrhythmogenesis. In the end-stages of disease, fibrosis replaces up to 50% of the myocardium and is a key determinant of patient outcome [9][10][11].
Although LVOT obstruction is common and easily treatable in HCM, a subset of patients without LVOT or LV cavity obstruction, referred to as nonobstructive HCM, develop progressive, symptomatic heart failure, despite guideline directed medical therapy, often leading to heart transplantation [12]. Histopathological features include myocyte hypertrophy, fibrosis, myocyte disarray and microvascular disease and are common with obstructive HCM. The genetic profiles of patients with nonobstructive HCM are also indistinguishable from those with obstructive HCM, raising questions about the sources of asymmetric hypertrophy. Human tissue for study of nonobstructive HCM can only be obtained postmortem or at the time of heart transplantation or other cardiac procedure.
Single cell and single nucleus transcriptomic analysis have facilitated an understanding of cellular phenotypes and interactions occurring in complex tissues such as the heart and high quality datasets of the normal human heart have been published [13][14][15]. We have recently performed single nucleus RNA-seq analysis of obstructive HCM and found significant alterations in intercellular communication pathways that affect the extracellular matrix, involving integrin-β1, thus providing a potential mechanism linking sarcomere dysfunction to extracellular matrix signaling [16]. We also compared the single nucleus RNA-seq profiles of obstructive and nonobstructive HCM and identified common and distinct pathogenic mechanisms [17]. Here, we analyze single nucleus transcriptomes in tissue from nonobstructive, end stage HCM in more detail and compare with normal tissue. As with obstructive HCM, we find that alterations in intercellular communication pathways affecting the extracellular matrix involve integrin-β1, but there are additional alterations in communication between fibroblasts, vascular cells, and dendritic cells not seen in obstructive HCM. In addition, there are increases in various molecular functions associated with adenylate cyclase signaling, ion channel function, protein phosphatase and protein kinase activity, suggesting a unique intercellular signaling milieu specific to nonobstructive HCM. These findings have important implications for the development of novel therapies for nonobstructive HCM.

Study Patients, Sample Collection and Processing
A total of 6 patients with end stage, nonobstructive HCM scheduled for cardiac transplantation were approached for written informed consent to allow their tissue to be used for research. Those who consented underwent surgery and tissue was collected. Explanted heart tissue sample processing was done as previously described [15]. 100 mg of collected interventricular septum tissue was minced into 1 mm 3 pieces, placed in 0.5 mL of CryoStor CS10 Freeze Media (STEMCELL Technologies), and stored in a MrFrosty (ThermoFisher) at 4°C for 10 minutes and then transferred to −80°C overnight. Bulk RNA was isolated from a piece of tissue using the Qiagen RNeasy Plus Micro kit and then assessed on the Agilent Bioanalyzer 2100. Samples with an RNA Integrity Number greater than 8.5 were used in library preparation. Sample collection was approved by the Tufts University/Medical Center Health Sciences Institutional Review Board under IRB protocol # 9487. Patient characteristics were obtained from the medical record and are shown in Table 1. Tissue from organ donor patients without underlying cardiac disease was obtained and processed as described previously [15,16].

Nuclei isolation, library preparation, and sequencing
Generation of single nuclei sequencing libraries was performed as previously described [15,16]. Cryopreserved samples were thawed at 37°C and placed on ice. Nuclei were isolated via Dounce homogenization as previously described [18]. Homogenates were filtered through a Pluristrainer 10 μM cell strainer (Fisher Scientific) into a pre-chilled tube. Nuclei were pelleted by centrifuging at 500 x g for 5 min at 4°C. Nuclei pellets were washed and pelleted according to manufacturer protocol (10x Genomics). Nuclei were stained with trypan blue and counted on a hemocytometer to determine concentration prior to loading of the 10x Chromium device and samples were diluted to capture ~10,000 nuclei. Nuclei were separated into Gel Bead Emulsion droplets using the 10x Chromium device according to the manufacturer protocol (10x Genomics). Sequencing libraries were prepared using the Chromium Single Cell 3´ reagent V2 kit according to manufacturer's protocol. Libraries were multiplexed and sequenced on a NovaSeq S4 (Illumina) to produce ~50,000 reads per nucleus. Single nuclei RNAseq data for normal IVS tissue from 4 donors is available in the Gene Expression Omnibus database under accession number GSE161921 [15]. Data for nonobstructive HCM IVS tissue and additional normal heart donors is available under accession number GSE181764.

Clustering of Cells by Gene Expression Pattern and Assignment of Cell Type Identity
Clustering of cells and assignment of cell identity was done as previously described [15][16][17]. Sequencing reads were processed using Cell Ranger version 6.0.1 [19]. The gene expression matrix was subset to only include reads from the nuclear genome. Quality control (QC) filtering, clustering, dimensionality reduction, visualization, and differential gene expression were performed using the R package Seurat version 3.0. Each dataset was filtered so that genes that were expressed in three nuclei or more were included in the final dataset. The dataset was further sublet to exclude nuclei that had fewer than 200 genes expressed to remove droplets containing only ambient RNA, and to exclude nuclei with greater than 2000 genes to remove droplets that contained two nuclei. Datasets were individually normalized and integrated using Seurat's SCTransform development workflow to reduce batch effects [20]. Optimal clustering resolution was determined using Clustree version 0.4.3 [21] to identify the resolution where the number of clusters stays stable and was determined to be 0.9 for the integrated dataset. Assignment of cell identity to each cluster was performed using four separate analyses. Expression of known cell-specific gene markers were used to identify major cell types, as done previously [13,[15][16][17]22]. The top 20-30 differentially expressed genes in each cluster were also compared with cell type gene expression markers from the PanglaoDB database https://panglaodb.se [23] to independently assign cell type. Entire sets of differentially expressed genes for each cluster were also subjected to Ingenuity Pathway Analysis [24] and their inferred functions were used to identify cluster cell types independently. Upregulated genes from each cluster were also subject to Gene Ontology biological process association using GoStats [25] and these associations were used to refine cell type assignment further.

Trajectory Analysis and Identification of Differentially Expressed Genes
Trajectory analysis was performed using Monocle3 [26] to determine the relationship between subtypes of cells identified in our clustering analysis as previously described [16,17]. Since our data do not represent a developmental time course, we determined the root nodes for each subtype by hierarchical clustering prior to generating trajectories and assigning pseudotime to each cell. Each cell type was analyzed in three cohorts: 1. Normal and nonobstructive HCM cells together; 2. Normal cells alone; 3. Nonobstructive HCM cells alone. Once trajectories were established, differential expression between Normal and HCM cells by cell type and cluster was performed by fitting a generalized linear regression model to each gene according to the formula log(yi) = β0 + βtxt. The coefficients β0 and βt were extracted from each model and tested for significant difference from zero using the Wald test, which assesses constraints on statistical parameters based on the weighted distance between the unrestricted estimate and its hypothesized value under the null hypothesis, where the weight is the precision of the estimate. A gene was determined to be differentially expressed between Normal and HCM conditions in a cell type or cluster if the Wald test produced an adjusted p-value ≤0.05.
Differentially expressed genes over trajectory paths in UMAP space (i.e., spatial autocorrelation) was performed in Monocle3 using Moran's I statistic. Moran's I statistic is a value that varies from −1 to 1, where −1 indicates perfect dispersion, 0 indicates no spatial autocorrelation, and 1 indicates perfect positive autocorrelation (i.e., nearby cells in have similar gene expression values in focal region of UMAP space). For each Normal and nonobstructive HCM cell type, a gene was determined to be differentially expressed over space if the associated Moran's I statistic value was positive, paired with a significant adjusted p-value ≤ 0.05, and expressed in ≥ 1% of associated cells. Since many genes showed differential expression over space, further conservative filtering was performed in which genes with Moran's I statistic available in a single class (i.e., Normal or HCM) were filtered by Moran's I statistic values >0.1. For genes with Moran's I statistics available in both classes (i.e., Normal and HCM), genes were filtered by an absolute difference >0.1. GO analysis of molecular function and biological process associated with differentially expressed genes was done using the online tools at uniprot.org/uniprotkb [27].

Analysis of Ligand-Receptor Pair Gene Expression to Discover Intercellular Communication Pathways
To quantify potential cardiac cell-cell communication in Normal and HCM hearts, cell communication networks were plotted in iGraph version 1.2.6 [28] and compared on the basis of ligand-receptor pair gene expression. Our cell-cell communication networks were derived as described previously [29], using a list of 2557 human ligand-receptor pairs [30] combined with another list of 3398 human ligand-receptor pairs [31], to give a total of 3627 unique human ligand-receptor pairs, largely as described previously [16,17]. A ligand or receptor for each cell type or cluster was considered expressed if the corresponding gene showed an above zero gene count in ≥ 20% of cells in our snRNAseq data, a threshold used commonly in other studies [16,17,29]. We initially analyzed the potential signaling interactions between the 8 cell types identified in our snRNAseq data. Lines in our cell networks connect two cell types and represent expressed human ligand-receptor pairs (i.e., potential cell-cell communication between a broadcasting (ligand) and recipient (receptor) cell types. Line color in our networks represents the broadcasting ligand source. Line thickness is proportional to the number of uniquely expressed ligand-receptor pairs. Cell-cell communication networks were also analyzed by fibroblast cluster along with other cell types, by cardiomyocyte cluster and other cell types and by fibroblast clusters and cardiomyocyte clusters. GO analysis of differentially expressed ligand receptor pairs was performed using the R package clusterProfiler [32].

Statistics
We used mixed effects models to analyze cell type specific differential expression, while taking into account the variability between and within subjects. A sample size of 6 cases and 6 controls, with an average of 3,000 cells per subject, will provide 80% power to detect fold change of expression ranging between 1.3 for an intracluster correlation f 0.01, and 2 for an intracluster correlation of 0.1. The power calculations were stable for sample sizes ranging between 500 and 3000 cells per subject. We used a Bonferroni correction for 10,000 tests to fix the level of significance. The power calculations were conducted using Power Analysis and Sample Size software (PASS).
Other statistical methods to cluster cells in UMAP space and control for batch effects using the built-in functionality of Seurat [20], to determine cluster stability by Clustree [21] and to perform Gene Ontology Analysis using well documented software programs [32] are described in above methods sections and in cited references. Statistical methods to compare gene expression along pseudotime trajectories using linear regression or spatial autocorrelation using the built-in functionality of Monocle3 are described above and in the original cited reference [26].

Study Approval
Explanted HCM heart sample collection was approved by the Tufts University/Medical Center Health Sciences Institutional Review Board under IRB protocol # 9487. All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki. Unused donor hearts were obtained in deidentified fashion from New England Donor Services under a Tufts University/Medical Center Health Sciences IRB approved research protocol after being designated as not human subjects research.

Patient characteristics
Explanted heart tissue was obtained from six HCM patients with severely symptomatic, nonobstructive, end stage disease who underwent cardiac transplantation. Patient characteristics are summarized in Table 1. Patients gave informed consent for their heart tissue to be used in research. The patients varied in age from 25 to 67. Five of six patients were female. Two of the patients carried pathogenic MYH7 mutations and the remaining four patients had no known mutations pathogenic for HCM. Four organ donors and datasets for the normal human IVS have been previously described [15] and are summarized in Table  2 with the two additional donors used in this study.

Nonobstructive HCM IVS tissue reveals extensive cardiomyocyte and fibroblast diversity
After sequencing and initial data processing with Cell Ranger software [19], each sample dataset was processed further to remove called nuclei that were likely droplets with ambient RNA, or droplets that contained two nuclei. The six HCM datasets and six donor heart datasets were combined into one dataset using the Seurat Integration function [20]. The final combined dataset included 49010 nuclei from nonobstructive HCM hearts and 34358 nuclei from donor hearts. Clustering of the integrated dataset revealed 22 cell populations within the IVS which was visualized using the dimensionality reduction algorithm uniform manifold approximation and projection (UMAP, Fig. 1A). Each point represents a single nucleus colored by cluster identity. Visualizing the integrated dataset by Normal and HCM datasets reveals that all 22 clusters are present in each condition and all samples contribute to all clusters ( Fig. 1B, C, D). No cell populations (i.e., clusters) were specific to either condition.
Cell identities were assigned to each cluster using known gene markers of expected cell types, differentially expressed genes queried against panglaoDB [23], gene ontology (GO) using GOstats [25], and Ingenuity Pathway Analysis [24]. Supplemental Table 1 lists the gene markers used to identify each cell type. Supplemental Table 2 lists the consensus cell identity assignment using all four methods. Upon assigning cell types, we found that similar cell types expressed similar markers and the associated clusters were expectedly positioned close to each other in UMAP space ( Fig. 2A, B). Interestingly, we see ten cardiomyocyte populations (i.e., clusters) and six fibroblast populations, revealing significant cardiomyocyte and fibroblast diversity in the IVS ( Fig. 2A, 2B), consistent with previous reports [15][16][17]. Other cell types identified included endothelial, smooth muscle, pericyte, neuronal, leukocyte, and dendritic cell populations. Among normal donor heart samples, assigned cardiomyocytes made up approximately 48.5% of the total cell population. Similarly, among nonobstructive HCM samples cardiomyocytes made up approximately 55.0% of the total cell population. Noncardiomyocytes make up the remaining cell population. Nonobstructive HCM does not appear to be associated with a shift in the relative proportion of cells in the heart.

Trajectory analysis and differential gene expression analysis reveal nonobstructive HCMassociated perturbations
Clustering analysis revealed diversity of both cardiomyocytes and fibroblasts. To gain insight into the potential relationships among the different cell populations for each cell type, clusters were projected onto pseudotime trajectory analysis using Monocle3 [26]. Trajectory analysis was performed on Normal cells only, nonobstructive HCM cells only, and both conditions combined for each of the 8 cell types identified in the above analysis. In single-cell trajectory analysis, a trajectory is a computed path that describes a cell type's biological progression through a dynamic process. Prior to building trajectory paths, the root node, or beginning, of each trajectory path was determined through hierarchical clustering of our single nuclei data by cell type. Then, trajectory paths for each assigned cell type (assigned in Seurat) were constructed in UMAP space using Monocle3 with Normal and HCM data together and individually. Once trajectory paths were established, pseudotime could be assigned to each cell. For all cell types, root nodes showed consistent placement among Normal only groups, HCM only groups, and Normal and HCM groups. Trajectory paths also showed similar basic structures among Normal only groups, HCM only groups, and Normal and HCM group (data not shown). Therefore, we were not able to distinguish any meaningful changes among trajectory paths. These findings suggest that the relationships between the various subtypes of each cell type do not vary significantly in nonobstructive HCM heart tissue in comparison to Normal heart tissue, as we have reported for obstructive HCM [16].
Establishment of trajectories facilitates the analysis of differential gene expression along the trajectory and between conditions along the trajectory. We analyzed differential gene expression for each cell type along their trajectories for the Normal and HCM populations. Differential expression analysis between Normal and HCM cells by assigned cell type and cluster was performed in which generalized linear regression models were fit for each gene. Resulting coefficients were tested for significant differences from zero, with an adjusted p-value cutoff ≤0.05 using the Wald Test. No differentially expressed genes were detected by this method.
Differential gene expression can also be analyzed along trajectory paths in UMAP space through spatial autocorrelation. Moran's I statistic and adjusted p-values were calculated for each gene in Normal only and HCM only groups within each cell type. When paired with a significant p-value (≤0.05), a Moran's I statistic value near zero indicated no spatial autocorrelation and a value near 1 indicated perfect positive autocorrelation in which a gene is expressed in a focal region of the UMAP space. Results showed between 115 and 6695 genes are differentially expressed along trajectory paths among cell types and between among cell types (Supplemental Table 3).
Since many genes showed differential expression over UMAP space, further conservative filtering was performed to identify genes of interest. Genes with Moran's I statistic available for a single condition, Normal or HCM, were filtered by a Moran's I statistic value above 0.1, while genes with Moran's I statistics available for both classes were filtered by an absolute difference >0.1. For each cell type, 1 to 18 genes passed the filter with 64 unique genes in total passing the filter among all cell types (Supplemental Table 4). The expression of these filtered genes was then plotted in UMAP space for their associated cell type in Normal and HCM conditions as shown for a subset of identified genes in cardiomyocytes (Supplemental Figure 1). Visual analysis of the UMAP plots revealed 41 genes with pronounced differences in their spatial expression between Normal and nonobstructive HCM conditions (summary in Table 3). Of these 41 genes, MYL2 is already known to be associated with human HCM. GO Enrichment analysis of these differentially expressed genes revealed significant enrichment for molecular functions involving actin binding, biological processes involving muscle development and muscle differentiation and cellular components including the sarcomere, myofibril and contractile fiber (Fig. 3A). Analysis of genes from this list showing increased expression in nonobstructive HCM revealed significant enrichment in molecular functions involving structural constituent of muscle and actin binding, no enriched biological processes and enrichment in cellular components involving the sarcomere, myofibril and contractile fiber (Fig. 3B). Analysis of genes from this list showing decreased expression in nonobstructive HCM revealed enrichment for molecular functions involving lipid transport and the ECM structural constituent, biological processes involving muscle differentiation, muscle development, lipid transport and lipid localization, and cellular components involving the endoplasmic reticulum membrane (Fig.  3C).

Ligand-receptor pair gene expression suggests alteration of intercellular communication in nonobstructive HCM
Nonobstructive HCM is associated with a general reduction in potential cell-cell communication but also increases in specific potential intercellular interactions-To quantify potential cardiac cell-cell communication in Normal and nonobstructive HCM tissue (n=405) was much lower than in Normal tissue (n=710, Fig.  4A, B). Analysis of the broadcast ligands by individual cell types demonstrates that 1) all cells except for smooth muscle cells show a broad decrease in both paracrine and autocrine ligand broadcasting in nonobstructive HCM (Fig. 4B, C, D) and 2) fibroblasts broadcast the greatest number of ligands to the other cell types in both Normal and nonobstructive HCM tissue. Analysis of receptor expression again showed a broad decrease across cell types, except for dendritic cells (Fig. 4B, C, D). Fibroblast communication is particularly high with fibroblast, endothelial, pericyte, smooth muscle, neuronal, and cardiomyocyte cells (Fig. 4). While intercellular communication was generally reduced in nonobstructive HCM tissue compared to Normal tissue, there was several cases of increased nonobstructive HCM communication, from fibroblasts, smooth muscle cells, pericytes and cardiomyocytes to dendritic cells and from smooth muscle cells to leukocytes (Fig 4C, D).
Changes in ligand-receptor pair gene expression imply alterations in molecular function-To assess the molecular functions potentially affected by changes in ligand-receptor signaling among our 8 identified cell types, GO enrichment analysis was performed on the ligands and receptors in expressed ligand-receptor pairs from both Normal and nonobstructive HCM tissue (Fig. 5). We observed a decrease in most identified ligand molecular functions (34 of 47) in the nonobstructive HCM heart tissue compared to the Normal heart tissue, with complete loss of functions involving insulin receptor binding, insulin-like growth factor receptor binding, lipoprotein particle receptor binding, and structural molecule activity conferring elasticity (note, these are low count normal MFs; Figure 5A). Large decreases in nonobstructive molecular functions relative to Normal (>100 count difference) were found for extracellular matrix structural constituent, extracellular matrix structural constituent conferring tensile strength (56%), growth factor binding (57%), peptidase regulator activity (75%), platelet-derived growth factor binding (56%), protease binding (74%; Figure 5A). Large increases in molecular functions in nonobstructive HCM tissue relative to normal tissue (>50 count difference) were noted for adenylate Codden

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript cyclase binding, calcium channel inhibitor activity, calcium channel regulator activity, channel inhibitor activity, protein kinase activator activity, protein N-terminus binding, protein phosphatase activator activity, protein serine threonine kinase activator activity and titin binding. When ligand molecular function is assessed in individual cell types, a general decrease in most nonobstructive molecular functions is observed compared to the normal condition ( Figure 5B-I) with generally more drastic decreases in cardiomyocytes, endothelial cells, dendritic cells, leukocytes, and neuronal cell ligand molecular functions ( Figure 5). Less drastic decreases and increases are observed in non-obstructive ligand molecular functions among fibroblasts, pericytes, and smooth muscle cells ( Figure 5).
Molecular functions associated with changes in receptor expression align with changes in ligand expression (Fig. 6A) including the finding of a general decrease in most nonobstructive molecular functions compared to the Normal condition (46 of 65 molecular functions; Figure 6). The largest decreases (>100 count difference) in nonobstructive molecular functions compared to Normal included amyloid beta binding, coreceptor activity, lipoprotein, peptide binding, and scavenger receptor activity ( Figure 6A). The largest increases (> 100 count difference) in nonobstructive molecular functions relative to Normal included amide binding, calcium molecular functions, calmodulin binding, channel/gated channel molecular functions, divalent inorganic cation transmembrane transporter activity, metal ion transmembrane transporter activity, protein kinase A activity, and sulfur compound binding ( Figure 6A). Trends in receptor molecular functions in individual cell types are notable for increases related to signal transduction in dendritic cells and leukocytes, particularly calcium channel activity and protein kinase A activity in dendritic cells.

Fibroblast subtypes show subtype-specific changes in potential cell-cell communication-To better understand the interaction between the different fibroblast
populations and other cell types in the IVS, we examined the potential ligand-receptor communication through cell network plots for all fibroblast subtypes alongside the other previously identified 7 cell types (Fig. 7A, C, D). In Normal tissue, all fibroblast subtypes broadcast ligands extensively to generate a dense intercellular network. In nonobstructive HCM tissue, the intercellular communication network for fibroblast subtypes shows a reduced number of broadcasting ligands compared to Normal tissue, (1502 versus 2636 respectively; Fig. 7A, B). Fibroblast cluster 2 generated the highest number of broadcasting ligands in both conditions. The largest decreases in potential interactions occurred between fibroblast cluster 2 and fibroblast clusters 2 through 6 (Supplemental Table 7) and between fibroblast cluster 5 and fibroblast clusters 2 through 6 (Supplemental Table 8). In fibroblast cluster 2 from nonobstructive HCM tissue, there was a disproportionate reduction in expression of cognate ligands for the ITGB1 receptor (COL18A1, COL5A1, COL5A2, FBLN1, FBN1, HSPG2, LAMC1, LGALS3BP, NID1, THBS2) and also of ligands for the LRP1 receptor (APP, C3, CALR, CTGF, HSPG2, LRPAP1, SERPINE2, SERPING1, TFPI). ITGB1 and LRP1 were expressed in both normal and nonobstructive HCM fibroblasts from cluster 2. Potential interaction between fibroblast cluster 2 and fibroblast clusters 3-5 was reduced in nonobstructive HCM primarily due to loss of the LRP1 receptor and several cognate ligands for the ITGB1 receptor (COL18A1, COL5A1, COL5A2, FBLN1, FBN1, HSPG2, LAMC1, LGALS3BP, NID1, THBS2). ITGB1 was expressed in fibroblast clusters

3-5 in both normal and nonobstructive HCM tissue.
Cardiomyocyte subtypes show subtype-specific changes in potential cell-cell communication-To better understand the potential interaction between the different cardiomyocyte populations and other cell types in the IVS, we performed L-R analysis for all ten cardiomyocyte subtypes and the other 7 cell types (Fig. 8). In nonobstructive HCM tissue, there is again a general decrease in potential intercellular communication (2634 L-R pairs in normal, reduced to 2026 in nonobstructive HCM; Fig. 8A, B), but the degree of reduction is lower than seen when subtypes are not considered or when fibroblast subtypes are considered. The largest decreases occur in ligands broadcast from endothelial cells and in receptors expressed by fibroblasts and the cardiomyocyte 4 cluster. Notably, there are increases in ligands broadcast by cardiomyocyte clusters 1, 2, 5, 9 and by smooth muscle cells and in receptors in cardiomyocyte clusters 1, 3 and 9 and dendritic cells in nonobstructive HCM (Fig. 8B). The greatest reductions in intercellular L-R connections involve fibroblasts broadcasting to fibroblasts and cardiomyocyte cluster 4 (Supplemental Table 9), and from endothelial cells broadcasting to fibroblasts, cardiomyocyte cluster 4 and cardiomyocyte cluster 8 (Supplemental Table 10). In both cases, reduction in intercellular communication is due to loss of LRP1 receptor expression in fibroblasts as well as loss of ligands for ITGB1 (COL11A1, FBLN1, FBN1, HSPG2, LAMB1, LAMC1, NID1, VCAN) in fibroblasts and endothelial cells. Notably, ITGB1 is expressed in fibroblasts and cardiomyocyte clusters 4 and 8 in both Normal tissue and nonobstructive HCM tissue. Increases in potential communication from cardiomyocyte clusters 2, 5, fibroblasts and smooth muscle cells to cardiomyocyte cluster 9 is due to gain of expression of ITGB1 in cardiomyocyte cluster 9 (Supplemental Table 11). Increased communication from cardiomyocyte cluster 5 to dendritic cells is also due to gain of expression of ITGB1 in dendritic cells (Supplemental Table 12).

Cardiomyocyte subtypes and fibroblast subtypes together show cell subtypespecific alterations in cell-cell communication-Given that cardiomyocytes and
fibroblasts are the most numerous cells in the heart and demonstrate the most subtypes in our study, we examined the intercellular networks between cardiomyocyte subtypes and fibroblast subtypes in normal tissue and nonobstructive HCM. The total number of expressed ligand-receptor pairs in the Normal tissue was greater than in the nonobstructive HCM tissue (3354 and 2399 expressed ligand-receptor pairs respectively, Fig. 9). The largest decreases in cellular communication involved fibroblast cluster 5, through reduction in both broadcasted ligands and expressed receptors. The largest reductions in L-R pairs occurred between fibroblast clusters 2 or 5 and fibroblast clusters 2-6 (Supplemental Table 13) and resulted from loss of ligands for ITGB1, which was still expressed in the recipient cells in both normal and nonobstructive HCM, or from loss of the LRP1 receptor in recipient fibroblast clusters except for fibroblast cluster 2. The largest increases in communication due to ligand broadcasting involved cardiomyocyte clusters 2, 5 and 9 and increases due to receptor expression involved cardiomyocyte clusters 1, 3, 9 and fibroblast cluster 1 (Fig.  9B). Increased communication to cardiomyocyte cluster 9 from fibroblast clusters 1-5 and cardiomyocyte clusters 2, 3 and 5 were due to gain of ITGB1 expression in cardiomyocyte cluster 9 (Supplemental Table 14).

Discussion
HCM has long been considered a disease of the sarcomere based on the association of sarcomere gene mutations with HCM, but the contributions from other loci are increasingly appreciated [2,3], as are potential mechanisms aside from sarcomere dysfunction [33]. We have previously examined single nuclei gene expression patterns in obstructive HCM. In this prior study, we did not identify any HCM-specific cell populations but observed a profound reduction in intercellular communication, especially in pathways mediated by ITGB1, and altered extracellular matrix gene expression, which may explain some of the nonmyocyte phenotypes seen in HCM such as fibrosis and link these phenotypes to alterations in mechanical signal transduction [16]. In the current study, we performed a single nucleus transcriptome analysis of tissue from patients with nonobstructive HCM, using tissue from the IVS harvested at the time of cardiac transplantation. This patient population is distinct from those with obstructive HCM in that they do not have LVOT obstruction and generally have more severe, intractable heart failure that results in cardiac transplantation. Determining gene expression patterns at the single cell level in nonobstructive HCM thus has the potential for facilitating precision medicine approaches for this specific population, through overlapping and unique pathogenic mechanisms relative to obstructive HCM.
The single nuclei gene expression patterns for nonobstructive HCM parallel those for obstructive HCM in that there are no disease-specific cell clusters and there is a reduction in intercellular communication, particularly in pathways involving ITGB1 and its associated ECM ligands. As with obstructive HCM, changes in gene expression involve sarcomere genes known to be associated with HCM and GO analysis indicates that gene expression changes are relevant to the extracellular matrix and sarcomere function, indicating that these are common, fundamental disease processes in HCM, regardless of LVOT obstruction. We have previously noted that a reduction in ECM associated gene expression and molecular function in obstructive HCM is counterintuitive based on the association between HCM and increased fibrosis and have suggested that the reduction in gene expression may represent a negative feedback loop that follows the accumulation of ECM protein in tissue [16]. Concordant findings in nonobstructive HCM can be explained similarly. Of note, a recent single cell RNA-seq study of enzymatically isolated cardiomyocytes from HCM patients undergoing myectomy surgery identified NPPA and XIRP2 as potential differentially regulated genes [34]. We did not observe differential expression of these genes in our current study or our previous studies [16,17], likely due to different methods. Precision medicine approaches that target integrin signaling and ECM interaction may one day prove useful for both obstructive and nonobstructive HCM.
A notable difference between obstructive and nonobstructive HCM in our current study, however, involves the unexpected increases in communication between fibroblasts, smooth muscle cells and pericytes with dendritic cells and between smooth muscle cells and leukocytes in the nonobstructive HCM condition relative to the control condition. Immune cells have been demonstrated to play important roles in development of cardiac hypertrophy and fibrosis in a pressure overload model of rodent hypertrophy (reviewed in [35]). Dendritic cells function as antigen presenting cells that activate T cell responses and depletion of these cells reduces cardiomyocyte hypertrophy, cardiac fibrosis, LV remodeling and inflammatory cell infiltration after pressure overload [36]. A role for dendritic cells in human heart failure is less well-established, although a trend toward elevated numbers in blood has been reported [37]. Integrin-mediated signaling has been shown to be an important determinant of dendritic cell function, with ITGB1 determining dendritic morphology [38] and promoting an anti-inflammatory phenotype in bone-marrow-derived dendritic cells [39]. Our data linking dendritic cells to fibroblasts, smooth muscle cells and pericytes by potential ITGB1 signaling is novel and suggests an important role for these cells in the pathological hypertrophy, cardiac fibrosis and vascular abnormalities seen in nonobstructive HCM. Our additional novel observation that signal transduction related molecular functions such as calcium channel activity and protein kinase A activity are increased in dendritic cells from nonobstructive HCM hearts also provides additional potential pathways for precisionmedicine guided therapy for nonobstructive HCM. These studies are consistent with a previously published but distinct analysis directly comparing obstructive and nonobstructive HCM [17].
We have previously noted extensive cardiomyocyte and fibroblast diversity in the human IVS, likely due to the diverse origins of the cells within this anatomic region [15][16][17]. We have also noted that the different cardiomyocyte and fibroblast subtypes exhibit specific intercellular communication profiles [16,17] and the phenomenon is observed again in this study. Specific cardiomyocyte and fibroblast subtypes appear to adopt distinct functional roles relevant to the disease process by shifting patterns of gene expression to increase or decrease interactions of cell types. An ongoing challenge is to determine which of these shifts is most relevant to the disease process and thus identify potential therapeutic targets. A limitation of our study, however, is that the relationship between single cell phenotypes and pathological anatomic features is unknown. Future spatial transcriptomic experiments assigning single cell populations to anatomic pathological features of HCM, such as areas of myocyte disarray, fibrosis, focal myocyte hypertrophy, and vascular abnormalities will likely be informative. Another limitation is that the alterations in ligand and receptor gene expression can infer but do not establish true intercellular interactions. Primary human tissue is not easily amenable to functional assessment of these ligand-receptor interactions. Future studies in model systems are needed to assess whether alterations in these ligandreceptor pathways affect the development of HCM. An additional challenge is the limited availability of tissue from patients with nonobstructive HCM and the small number of patients in our study, which precludes comparison of sarcomere gene mutation positive and negative patients and limits generalizability. Nevertheless, our study is among the first to characterize transcription patterns in nonobstructive HCM patients at single cell resolution and provides an important resource for future investigation of the complex cellular interplay in nonobstructive HCM.     A. GO enrichment analysis of all differentially expressed genes in HCM from Table 3. B. GO enrichment analysis of Table 3 genes increased in nonobstructive HCM. C. GO enrichment analysis of Table 3 genes decreased in nonobstructive HCM.