Single-cell RNA-sequencing Unravels Heterogeneity of Cardiomyocytes and Signaling Pathways in Pressure Overload

doi:10.21203/rs.3.rs-37034/v1

Download PDF

Research

Single-cell RNA-sequencing Unravels Heterogeneity of Cardiomyocytes and Signaling Pathways in Pressure Overload

https://doi.org/10.21203/rs.3.rs-37034/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

This study aimed to unravel the heterogeneity of cardiomyocytes and probed out hub genes and hub pathways for cardiac hypertrophy based on transverse aortic constriction (TAC) mouse models using single-cell RNA sequencing (scRNA-seq).

Methods

scRNA-seq data of TAC mouse models were retrieved from the GSE95140 dataset. After filtering, cell clusters were detected using scRNA-seq data, followed by identification of differentially expressed genes (DEGs). Then, functional enrichment analysis of DEGs was presented. GSVA scores of hub pathways were calculated. After that, hub genes were detected by protein-protein interaction (PPI) network and expression association analysis. Cell subtypes were clustered using UMAP and the expression patterns of hub genes across different cell subtypes and different stages of cardiac hypertrophy were visualized. Finally, hub genes and hub pathways were verified using the GSE76 and GSE36074 datasets.

Results

Following data filtering and normalization, 3408 DEGs were identified between TAC and sham operation. As shown functional enrichment analysis, hub pathways were identified including cardiac hypertrophy, ion transport, myocardial remodeling, apoptosis, HIF pathway and metabolise. Eight hub genes (Vldlr, Ugp2, Tgm2, Pygm, Flnc, Ctsd, Clu and Atp1b1) with the highest degree in the PPI network and the strongest correlation with GSVA calculated score of hub pathways were identified for cardiac hypertrophy. Six cell subtypes were clustered, composed of fibroblast, CM-A, CM-V, trabecular CM and endothelial cell. There was a distinct heterogeneity in the expression patterns of hub genes and the GSVA scores of hub pathways across different cell clusters and different stages of cardiac hypertrophy. The hub genes and hub pathways were externally verified by the two independent datasets.

Conclusion

Our findings identified hub genes and hub pathways for cardiac hypertrophy, which had a distinct heterogeneity across different cell clusters and different stages of cardiac hypertrophy.

Translational Medicine

Cardiac hypertrophy

Transverse aortic constriction

Single-cell RNA-sequencing

Signaling Pathways

Persistent cardiac hypertrophy is the dominating cause of heart failure [1]. Generally, it has been believed that cardiac hypertrophy can be a crucial therapeutic target for chronic heart failure [2]. Cardiac hypertrophy is usually characterized by an increase in the size of myocardial cells and thickening of the ventricular wall. Continuous hypertrophic stimulation (such as pressure overload [3], ischemia [4], and hypoxia [5]) gradually transforms the compensatory response into irreversible pathological myocardial hypertrophy, leading to cardiomyocyte apoptosis, fibrosis of extracellular matrix and abnormal expression of cardiac fetal gene expression [6]. At present, there is no effective drug therapy that can reverse the progression from pathological cardiac hypertrophy to heart failure. It is widely accepted that active intervention in the early stages of heart failure is an effective method to suppress pathological cardiac hypertrophy [7]. Despite continuous improvements in the diagnosis and treatment of cardiac hypertrophy, the mortality is close to 25% to 50% within 5 years following diagnosis [8]. Elucidation concerning the complex signaling mechanisms during myocardial hypertrophy may accelerate the improvement in treatment, thereby improving the quality of life of patients with cardiac hypertrophy [9]. Thus, it is urgent and necessary to probe out novel and effective therapeutic targets for preventing and cutting down pathological cardiac hypertrophy.

A single cardiomyocyte is the basic unit of gene regulation. Gene expression is the basis for determining the phenotype of cardiomyocytes and cardiac function, but it is still unclear which genetic programs are involved in maintaining and destroying the steady state of the heart. Compared to large amounts of RNA-seq that only provide an average expression signal of millions of cells, scRNA-seq can simultaneously analyze more than 10,000 single-cell transcriptomes, thereby characterizing novel cell clusters. Thus, scRNA-seq can reliably identify individual cell subpopulations and reveal the unique changes foreach cell type. Also, it can elucidate the heterogeneity of gene expression patterns in cardiomyocyte populations between healthy and cardiac hypertrophy.

Transverse aortic constriction (TAC) surgery directly reduces the inner diameter of the aortic arch and increases the left ventricular afterload, which is widely utilized to research pathological cardiac hypertrophy as well as heart failure [10-12]. In this study, we identified individual cardiomyocyte subtypes and probed out hub genes and hub pathways for cardiac hypertrophy by TAC scRNA-seq. Furthermore, we elucidated the heterogeneity of gene expression patterns at different stages of cardiac hypertrophy. Thus, our findings offered novel insights into the mechanisms of cardiac hypertrophy and underlying therapeutic targets.

scRNA-seq acquisition and preprocessing

Single-cardiomyocyte RNA-seq data were retrieved from GEO database (accession: GSE95140 dataset; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE95140) [13]. There were 396 single-cardiomyocyte transcriptomes from mice at the third (D3) day, the first week (W1), the second week (W2), the fourth week (W4) and the eighth week (W8) following TAC or sham operation in the GSE95140 dataset on the platform of GPL17021. Furthermore, microarray data in the GSE76 dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE76) were also obtained from GEO database based on the GPL32 platform, composed of mice at 1 hour, 4 hour, 24 hour, 48 hour, 1 week and 8 week after pressure-overload induced cardiac hypertrophy or sham operation. Microarray expression profile in the GSE36074 dataset downloaded from GEO database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36074) included seven mice with only hypertrophy, seven mice that had evident heart failure following pressure overload and seven mice treated with sham operation on the GPL1261 platform [14].

scRNA-seq clustering by Seurat

The “DropletUtils” package was utilized to detect the expression of each cell. Cells without any expression were filtered out. By the calculating “QC-Metrics” function in the “scater” package, the gene expression in the cell was counted [15]. Under the screening conditions of mitochondrial gene≤5% and ribosomal gene≥10%, the cells was further filtered. The expression matrix of each sample after filtering were normalized by “NormalizeData” function in the Seurat package (version 3.0) [16]. Through the “FindVariableFeatures” function in the Seurat package, the genes with the most obvious differences among cells were selected. The linearly scale of the expression data was presented using the “ScaleData” function in the Seurat package. Then, principal component analysis (PCA) was carried out by the “RunPCA” function in the Seurat package. Principal components (PCs) with >70% large standard deviations were selected. Using “RunUMAP” in the Seurat package, UMAP dimensionality reduction analysis was performed. By the “FindAllMarkers” function in the “Seurat” package, differentially expressed genes between different cells, with the threshold of log |fold change (FC)|>0.5, the expression ratio of the cell population>0.25, and P≤0.05), thereby obtaining marker genes.

Gene Set Enrichment Analysis (GSEA)

GSEA was presented to evaluate microarray expression data at a group of gene sets [17]. The enrichment score (ES) was firstly calculated utilizing the Kolmogorov–Smirnov algorithm. Then, we estimated the statistical significance of ES in line with the empirical phenotype replacement test procedure. The significance level was adjusted by adjusted multiple hypothesis tests The ES for each gene set was normalized into the enrichment score (NES). The false discovery rate (FDR) corresponding to each NES was calculated.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis

GO enrichment analysis was presented using “Clusterprofiler” package in R [18]. KEGG pathway map was visualized according to Metascape online database, a gene annotation and analysis resource (http://metascape.org) with the threshold of overlap≥3, P-value≤0.05 and enrichment score<1.5 [19].

Gene Set Variation Analysis (GSVA)

The activation of gene sets was analyzed by GSVA utilizing the “GSVA” package in R. GSVA is a non-parametric and unsupervised algorithm for estimating the enrichment score of hub pathways based on each cell transcription expression matrix in the TAC model. In each cell, after assigning different groups, GSVA scores were calculated by the genefu package in R, which were compared by the one-way ANOVA between multiple groups.

Protein-protein interaction (PPI) network

Genes enriched in hub pathways were extracted. Using STRING online database (version 11.0; https://string-db.org/), a PPI network was conducted with the threshold of 0.4 [20]. The degree of each node was calculated utilizing Cytoscape software (version 3.8.0; https://cytoscape.org/) [21]. Spearson correlation analysis between GSVA calculated scores of hub pathways and expression levels of enriched genes was presented. Correlation coefficient>0.4 as well as P<0.05 was considered significantly correlated. For each hub pathway, a gene with the highest degree in the PPI network and the strongest correlation with GSVA calculated score of hub pathway was screened out as a hub gene.

External verification of hub genes and hub pathways

The expression levels of hub genes and GSVA calculated score of hub pathways were validated using two independent datasets including GSE76 and GSE36074 datasets.

scRNA-seq clustering by Seurat for cardiac hypertrophy

An overview of flowchart in this study is shown in Figure 1A. Single-cardiomyocyte RNA-seq data were obtained from GSE95140 dataset. GSEA was used to analyze cardiac hypertrophy-related KEGG pathways. As shown in Figure 1B, citrate cycle TCA cycle (ES=0.84, NES=1.82, p-value<0.01, Size=30), cysteine and methionine metabolism(ES=0.60, NES=1.78, p-value<0.01, Size=32), fatty acid metabolism(ES=0.85, NES=2.21, p-value<0.01, Size=39), lysine degradation(ES=0.73, NES=2.23, p-value<0.01, Size=41), propanoate metabolism(ES=0.81, NES=2.08, p-value<0.01, Size=31), apoptosis (ES=-0.33, NES=-1.20, p-value=0.04, Size=80), arrhythmogenic right ventricular cardiomyopathy (ARVC; ES=-0.37, NES=-1.34, p-value=0.03, Size=74), dilated cardiomyopathy (ES=-0.35, NES=-1.30, p-value=0.05, Size=85), ECM receptor interaction (ES=-0.35, NES=-1.30, p-value=0.04, Size=82) and focal adhesion (ES=-0.40, NES=-1.50, p-value<0.01, Size=192) were significantly enriched by differentially expressed genes between TAC and sham groups. Furthermore, cardiac muscle contraction (ES=0.65, p-value<0.01, Size=69), citrate cycle TCA cycle (ES=0.79, p-value<0.01, Size=30), fatty acid metabolism (ES=0.70, p-value<0.01, Size=32), glutathione metabolism (ES=0.68, p-value<0.01, Size=22) and glycolysis gluconeogenesis (ES=0.76, p-value<0.01, Size=46) could be significantly related to the development of cardiac hypertrophy according to TAC time series.

Cardiomyocytes were screened for further analysis (Figure 2A; Table S1). Then, PCA was used to select PC with large standard deviations. Elbow diagram showed that PC was set as 12 (Figure 2B). UMAP was used to classify these cardiomyocytes into six clusters (Figure 2C). After normalization, 3408 highly variable genes were screened and the top 20 genes were visualized, as shown in Figure 2D. Heat maps depicted the top ten marker genes for each cluster (Figure 2E).

Identification of hub pathways for cardiac hypertrophy

Consequently, a total of 288 markers were detected based on Wilcoxon signed rank test among the different cell clusters (Table S2). To probe out potential biological functions of differentially expressed genes, we presented functional enrichment analyses. GO enrichment analysis results showed that these genes were primarily enriched in cardiac hypertrophy-related biological processes including muscle cell development (Count=22, p-value=1.57E-14), myofibril assembly (Count =15, p-value=4.58E-14), ribonucleotide metabolic process(Count =28, p-value=5.57E-14), ribose phosphate metabolic process (Count =28, p-value=9.14E-14) and cardiac muscle tissue development (Count =22, p-value=2.32E-13). Furthermore, these genes participate in myocardial cell components such as contractile fiber (Count =25, p-value=7.23E-18), myelin sheath (Count =23, p-value=2.14E-14), mitochondrial inner membrane (Count =27, p-value=8.43E-14) and actin cytoskeleton (Count =20, p-value=3.95E-07). Also, they had several crucial molecular functions like ubiquitin protein ligase binding(Count =15, p-value=4.20E-05), actin binding(Count =16, p-value=1.86E-04), ion channel binding(Count =9, p-value=3.02E-04), iron ion binding (Count =5, p-value=3.02E-04) and antioxidant activity(Count =7, p-value=6.84E-04) (Figure 3A; Table S3). As visualized in the KEGG pathway map, these genes were mainly involved in cardiac hypertrophy, ion transport, myocardial remodeling, apoptosis, HIF pathway and Metabolize (Figure 3B; Table S4). As shown in Figure 3C, there were distinct differences in the GSVA calculated scores of these KEGG pathways at different stages of cardiac hypertrophy compared to sham operation. Concretely, compared to sham group, cardiac hypertrophy scores were significantly higher in TAC D3, W1, W2 and W4 groups (all p-value<0.05). However, no statistical difference was detected between sham and TAC W8 groups. Compared to sham group, there were significantly higher apoptosis scores in TAC D3, W1and W4 groups (all p-value<0.05) not TAC W2 and W8 groups. Lower metabolise scores were found in TAC W1 and W8 groups (both p-value<0.05) not TAC D3, W2 and W4 groups than sham group. Also, our results showed that HIF pathway scores were distinctly lower in TAC W1 and W2 groups (both p-value<0.05) not TAC D3, W4 and W8 groups than sham group. Furthermore, compared to sham group, there were distinctly lower ion transport scores in TAC W1 group (p-value<0.05), not TAC D3, W2, W4 and W8 groups. Myocardial remodeling scores were remarkably higher in TAC D3, W1, W2, W4 and W8 groups than sham group (all p-value<0.05). Figure 3D depicted that the difference of the GSVA calculated scores of these KEGG pathways across different cell clusters.

Identification of hub genes and their expression in each cell subtype for cardiac hypertrophy

Eight hub genes with the highest degree in the PPI network and the strongest correlation with GSVA calculated score of hub pathways were identified for cardiac hypertrophy, including Ugp2 in the metabolise pathway (Correlation coefficient=0.47, Degree=35), Atp1b1 in the ion transport pathway(Correlation coefficient=-0.43, Degree=12), Flnc and Tgm2 in the cardiac hypertrophy pathway (Flnc: Correlation coefficient=0.48, Degree=22; Tgm2: Correlation coefficient=0.43, Degree=17), Ctsd and Clu in the apoptosis pathway (Ctsd: Correlation coefficient=0.45, Degree=20; Clu: Correlation coefficient=0.42, Degree=18), Pygm and Vldlr in the HIF pathway (Pygm: Correlation coefficient=0.47, Degree=16 Vldlr: Correlation coefficient=0.50, Degree=14) and Flnc in the myocardial remodeling (Correlation coefficient=0.46, Degree=17), as shown in Figure 4A and Table S5. Then, the six cell clusters were defined as fibroblast, CM-A, CM-V, trabecular CM and endothelial cell by singleR package (Figure 4B).

Additionally, Figure 4C showed that there was a high heterogeneity in expression of the eight hub genes across different cell clusters. For example, Vldlr was up-regulated in CM-A (average expression>1.5 and percent expressed=100) and down-regulated in trabecular CM and CM-V (average expression<-1.0 and percent expressed=100). Compared to trabecular CM and CM-A, there was a low percentage of Tgm2 expression in fibroblast, endothelial cell and CM-V (average expression>1 and percent expressed<99.0%).

Furthermore, we compared the expression patterns of these hub genes at different stages of cardiac hypertrophy in mice. Compared to sham operation, the expression levels of Atp1b1 and Vldlr were significantly elevated in mice treated by TAC at the D3, W1, W2 and W4 (all p-value<0.05; Figure 4D). For Clu, there were distinctly higher expression levels in mice treated by TAC at the D3, W1, W4 and W8 than sham operation group (all p-value<0.05). Ctsd expression had an increase expression levels in mice treated by TAC at the W2, W4 and W8 (all p-value<0.05). Compared to sham operation, both Flnc and Tgm2 expression was notably elevated in mice with TAC at the D3, W1, W2 and W8 (all p-value<0.05). Moreover, the expression levels of Pygm and Ugp2 were prominently higher in mice with TAC at the D3, W1, W2, W4 and W8 than in mice with sham operation (all p-value<0.05). These results indicated that there was a heterogeneity of expression patterns hub genes at different stages of cardiac hypertrophy.

Validation of hub gene expression and hub pathways using external datasets

Hub genes and hub pathways were validated using two external datasets (GSE76 and GSE36074). In the GSE76 dataset, our results showed that the expression levels of Atp1b1, Clu, Ctsd, Flnc, Pygm, Tgm2 and Ugp2 were distinctly higher in mice at 1w and 8w following pressure-overload induced cardiac hypertrophy than those in mice with sham operation for 1w (all p-value<0.05; Figure 5A). In converse, Vldlr had evidently lower expression levels in mice at 1w and 8w following pressure-overload induced cardiac hypertrophy than those in mice with sham operation for 1w (both p-value<0.05). As expected, there was no significant difference in their expression between mice with sham operation for 1w and for 8w. Similar to the validation results using the GSE76 dataset, the expression of Atp1b1, Ctsd, Flnc, Pygm, Tgm2 and Ugp2 were significantly increased and the expression of Vldlr was decreased both in mice with heart failure and only hypertrophy induced by TAC from the GSE36074 dataset (all p-value<0.05). For Clu, there was higher expression level in mice with heart failure (p-value<0.05) not only hypertrophy than mice with sham operation (Figure 5B).

The GSVA calculated scores of hub pathways were calculated based on the GSE76 and GSE36074 datasets. As shown in Figure 5C, the GSVA calculated scores of apoptosis, ion transport, metabolise and myocardial remodeling were obviously higher in mice at the first and eighth week after induced by pressure overload compared to sham operation group (all p-value<0.05). Cardiac hypertrophy had a higher GSVA calculated score in mice at the first week after induced by pressure overload than in mice with sham operation for 1w (p-value<0.05). For HIF pathway, its lower GSVA calculated score was found in mice at the first week after induced by pressure overload compared to mice with sham operation for 1w (p-value<0.05). In the GSE36074 dataset, the GSVA calculated scores of apoptosis, cardiac hypertrophy, metabolise and myocardial remodeling were prominently higher both in mice with heart failure and only hypertrophy induced by TAC (all p-value<0.05; Figure 5D). Conversely, the GSVA calculated score of HIF pathway was significantly decreased both in mice with heart failure and only hypertrophy induced by TAC (both p-value<0.05). Furthermore, ion transport had a distinctly lower GSVA calculated score in mice with heart failure compared to sham operation (p-value<0.05).

In this study, we identified hub pathways and hub genes for cardiac hypertrophy using scRNA-seq, which were verified by two independent datasets. Furthermore, six cell subtypes were clustered, composed of fibroblast, CM-A, CM-V, trabecular CM and endothelial cell. A distinct heterogeneity in the expression patterns of hub genes and the GSVA scores of hub pathways was found across different cell clusters and different stages of cardiac hypertrophy. These findings revealed novel insights into the molecular mechanisms of cardiac hypertrophy and these hub genes could become potential therapeutic targets for cardiac hypertrophy.

Increasing evidence suggests that it is of importance to probe out the mechanisms of cardiac hypertrophy by using scRNA-seq analysis. For instance, a previous study found that the abnormal expression of Myh7 in cardiomyocytes following pressure overload was regulated by time through scRNA-seq [22]. Another study reported that cardiac differentiation of human pluripotent stem cells could depend on HOPX based on scRNA-seq analysis [23]. Furthermore, Michail Yekelchyk et al revealed the heterogeneity between cardiomyocytes in stress-induced myocardial hypertrophy utilizing scRNA-seq [24]. In this study, we identified eight hub genes (including Vldlr, Ugp2, Tgm2, Pygm, Flnc, Ctsd, Clu and Atp1b1) with the highest degree in the PPI network and the strongest correlation with GSVA calculated score of hub pathways for cardiac hypertrophy. Na, K-ATPase, also known as sodium pump, contains α- and β-subunits, which served as an ion pump and a signal transducer [25]. The enzyme activity of Na, K-ATPase and its α1, α 3 and β 1 subtypes have been confirmed to be decreased in the human heart failure. In this study, we found that the expression of Atp1b1 (as called Na, K-ATPase β 1) was distinctly elevated in mice at 1w and 8w following pressure-overload induced cardiac hypertrophy compared to sham operation mice. Moreover, its high expression was detected in mice with heart failure or only hypertrophy induced by TAC compared to sham operation. In mice with Atp1b1 knockdown, cardiac hypertrophy as well as myocardial contractility disorder occurred after modeling TAC [26]. Clu (as called Clusterin), has been considered as a biomarker for heart failure. Recently, Clu expression could be activated by IGF1-PI3K pathway both in heart tissues and cardiomyocytes, which is involved in physiological cardiac hypertrophy as well as myocardial protection [27]. We found that Clu expression was higher in mice at the first and eighth week after induced by pressure overload. Consistent with previous study, in mice with heart failure, the expression of Clu was notably elevated compared to sham operation. However, there was no significant difference between mice with only cardiac hypertrophy and sham operation. Our findings also indicated that Clu possesses potential value as a marker for cardiac hypertrophy. Ctsd (alias: Cathepsin D) is a lysosomal protease that maintains the homeostasis of cardiomyocytes [28]. Ctsd has been found to be down-regulated both in myocardial tissue and blood of patient with cardiac hypertrophy [29]. Conversely, our results showed that Ctsd had remarkably higher expression levels in mice induced by pressure overload after 1 week and 8 weeks. Both in mice with heart failure and with only cardiac hypertrophy, the expression of Ctsd was prominently increased. Flnc (as known as Filamin C), an actin cross-linked protein, is specifically expressed in myocardial tissue, which is crucial to maintain the integrity of myocardial structure [30]. Positive expression of Flnc could response to pressure overload-induced cardiac hypertrophy [31], which was consistent with our findings. Pygm (Glycogen Phosphorylase, Muscle Associated) expression was distinctly higher in mice induced by pressure overload following 1 week and 8 weeks. Furthermore, both in mice with heart failure and with only cardiac hypertrophy, Pygm expression was evidently elevated. Previously, Tgm2 (Transglutaminase 2) expression is upregulated in myocardial tissues with pressure overload, which participates in promoting ventricular diastolic dysfunction and preventing ventricular dilation [32]. Activation of Tgm2 could contribute to age-related diastolic dysfunction [33]. Similarly, our results suggested that Tgm2 expression was distinctly activated in mice at the first and eighth week after pressure overload. Additionally, its expression was elevated in mice with heart failure or with only cardiac hypertrophy. There is no research concerning about Ugp2 (UDP-Glucose Pyrophosphorylase 2) on cardiac hypertrophy. In this study, Ugp2 had significantly higher expression levels both in mice at the first and eighth week after pressure overload. In addition, it was validated that Ugp2 expression was elevated in mice with heart failure or with only cardiac hypertrophy. In converse, Vldlr (very low-density lipoprotein receptor) expression was reduced both in mice at the first and eighth week after pressure overload compared to sham operation. Also, its low expression was detected in mice with heart failure or with only cardiac hypertrophy.

In this study, we identified six hub pathways for cardiac hypertrophy, including cardiac hypertrophy, ion transport, myocardial remodeling, apoptosis, HIF pathway and metabolise. Intriguingly, we found that the activity of apoptosis, cardiac hypertrophy, ion transport, metabolise and cardiac remodeling was distinctly increased in mice induced by pressure overload following 1 week. Conversely, HIF pathway activity was reduced in mice at the first week after pressure overload. Pressure overload could promote apoptosis of cardiomyocytes [34]. Several drugs (such as Tamarixetin [35], growth hormone releasing hormone [36] and Resveratrol [5]) have been found to inhibit cardiomyocyte apoptosis by targeting apoptosis-related genes. Increasing evidence suggests that cardiac ion channels such as voltage-gated Ca²⁺, Na⁺, K⁺ channels are regulated by multiple factors and their loss of control could induce arrhythmias and myocardial systolic disorders [37]. Metabolic changes play a key role in cardiac remodeling. Metabolic-mediated gene expression changes and signal transduction are essential to the physiological growth of the heart, while metabolic inefficiency and disorders of coordinated anabolic activity can induce cardiac remodeling cardiac remodeling [38]. Unlike persistent hypoxia exposure or strenuous exercise, severe hypoxia may place a premium on cardiomyocyte toxicity and hypertrophy [39]. It has been reported that Apigenin can reduce myocardial hypertrophy caused by pressure overload by down-regulating HIF-1α in rat heart [40].

Our research has some limitations. scRNA-seq cannot comprehensively evaluate all transcripts in cells, but is limited to 15%~30% of transcriptome. Moreover, the hub genes and hub pathways identified by our study are still lack of clinical and basic experimental verification. However, our research clearly suggests that there is a distinct heterogeneity between cardiomyocytes induced by pressure overload. Additionally, these hub genes and hub pathways were verified by two independent datasets to ensure the accuracy of the results. Therefore, our research will provide valuable resources for future research on the mechanism of cardiac hypertrophy.

In this study, we comprehensively analyzed scRNA-seq data from mice induced by TAC at the D3, W1, W2, W4 and W8. Six cell subtypes were clustered, composed of fibroblast, CM-A, CM-V, trabecular CM and endothelial cell. Furthermore, we identified hub pathways (including cardiac hypertrophy, ion transport, myocardial remodeling, apoptosis, HIF pathway and metabolise) and hub genes (including Vldlr, Ugp2, Tgm2, Pygm, Flnc, Ctsd, Clu and Atp1b1) for cardiac hypertrophy, which had a distinct heterogeneity across different cell clusters and different stages of cardiac hypertrophy. Thus, these hub genes and hub pathways deserve more in-depth research.

TAC: transverse aortic constriction; scRNA-seq: single-cell RNA sequencing (scRNA-seq; DEGs: differentially expressed genes; PPI: protein-protein interaction; PCA: principal component analysis; PCs: Principal components; GSEA: Gene Set Enrichment Analysis; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; GSVA: Gene Set Variation Analysis.

Availability of data and materials

All of the data were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database, and a research ethics application was not needed for this study.

Ethics approval and consent to participate

Not applicable

Competing interests

The authors have no conflict of interest to disclose.

Consent for publication

Not applicable.

Funding

This work was supported by the Guangzhou high-level Clinical Pivot Subjects Construction project, China (011102-02) and grant from Key-Area Research and Development Program of Guangdong Province(2019B020227001).

Acknowledgements

Not applicable.

Authors’ contributions

RZ, WS and MZ: takes responsibility for all aspects of the reliability and freedom from bias of the data presented and their discussed interpretation, drafting the article. WC, WL, NZ and LP: takes responsibility for statistical analyses, and interpretation of data. LM and XC: takes responsibility for full text evaluation and guidance, final approval of the version to be submitted. All authors read and approved the final manuscript.

Chen C, Zou LX, Lin QY, Yan X, Bi HL, Xie X, Wang S, Wang QS, Zhang YL, Li HH: Resveratrol as a new inhibitor of immunoproteasome prevents PTEN degradation and attenuates cardiac hypertrophy after pressure overload. Redox Biol 2019, 20:390-401.
Yu Q, Kou W, Xu X, Zhou S, Luan P, Xu X, Li H, Zhuang J, Wang J, Zhao Y et al: FNDC5/Irisin inhibits pathological cardiac hypertrophy. Clin Sci (Lond) 2019, 133(5):611-627.
Wang C, Yuan Y, Wu J, Zhao Y, Gao X, Chen Y, Sun C, Xiao L, Zheng P, Hu P et al: Plin5 deficiency exacerbates pressure overload-induced cardiac hypertrophy and heart failure by enhancing myocardial fatty acid oxidation and oxidative stress. Free Radic Biol Med 2019, 141:372-382.
Omidkhoda N, Wallace Hayes A, Reiter RJ, Karimi G: The role of MicroRNAs on endoplasmic reticulum stress in myocardial ischemia and cardiac hypertrophy. Pharmacol Res 2019, 150:104516.
Guan P, Sun ZM, Wang N, Zhou J, Luo LF, Zhao YS, Ji ES: Resveratrol prevents chronic intermittent hypoxia-induced cardiac hypertrophy by targeting the PI3K/AKT/mTOR pathway. Life Sci 2019, 233:116748.
Schauer A, Adams V, Poitz DM, Barthel P, Joachim D, Friedrich J, Linke A, Augstein A: Loss of Sox9 in cardiomyocytes delays the onset of cardiac hypertrophy and fibrosis. Int J Cardiol 2019, 282:68-75.
Lu X, He Y, Tang C, Wang X, Que L, Zhu G, Liu L, Ha T, Chen Q, Li C et al: Triad3A attenuates pathological cardiac hypertrophy involving the augmentation of ubiquitination-mediated degradation of TLR4 and TLR9. Basic Res Cardiol 2020, 115(2):19.
Zhao Y, Wang C, Hong X, Miao J, Liao Y, Hou FF, Zhou L, Liu Y: Wnt/β-catenin signaling mediates both heart and kidney injury in type 2 cardiorenal syndrome. Kidney Int 2019, 95(4):815-829.
Oldfield CJ, Duhamel TA, Dhalla NS: Mechanisms for the transition from physiological to pathological cardiac hypertrophy. Can J Physiol Pharmacol 2020, 98(2):74-84.
Byun J, Del Re DP, Zhai P, Ikeda S, Shirakabe A, Mizushima W, Miyamoto S, Brown JH, Sadoshima J: Yes-associated protein (YAP) mediates adaptive cardiac hypertrophy in response to pressure overload. J Biol Chem 2019, 294(10):3603-3617.
Shang L, Weng X, Wang D, Yue W, Mernaugh R, Amarnath V, Weir EK, Dudley SC, Xu Y, Hou M et al: Isolevuglandin scavenger attenuates pressure overload-induced cardiac oxidative stress, cardiac hypertrophy, heart failure and lung remodeling. Free Radic Biol Med 2019, 141:291-298.
Zhang B, Zhang P, Tan Y, Feng P, Zhang Z, Liang H, Duan W, Jin Z, Wang X, Liu J et al: C1q-TNF-related protein-3 attenuates pressure overload-induced cardiac hypertrophy by suppressing the p38/CREB pathway and p38-induced ER stress. Cell Death Dis 2019, 10(7):520.
Nomura S, Satoh M, Fujita T, Higo T, Sumida T, Ko T, Yamaguchi T, Tobita T, Naito AT, Ito M et al: Cardiomyocyte gene programs encoding morphological and functional signatures in cardiac hypertrophy and failure. Nat Commun 2018, 9(1):4435.
Skrbic B, Bjørnstad JL, Marstein HS, Carlson CR, Sjaastad I, Nygård S, Bjørnstad S, Christensen G, Tønnessen T: Differential regulation of extracellular matrix constituents in myocardial remodeling with and without heart failure following pressure overload. Matrix Biol 2013, 32(2):133-142.
McCarthy DJ, Campbell KR, Lun AT, Wills QF: Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R. Bioinformatics 2017, 33(8):1179-1186.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R: Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol 2018, 36(5):411-420.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES et al: Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005, 102(43):15545-15550.
Yu G, Wang LG, Han Y, He QY: clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 2012, 16(5):284-287.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK: Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun 2019, 10(1):1523.
von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B: STRING: a database of predicted functional associations between proteins. Nucleic Acids Res 2003, 31(1):258-261.
Doncheva NT, Morris JH, Gorodkin J, Jensen LJ: Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res 2019, 18(2):623-632.
Satoh M, Nomura S, Harada M, Yamaguchi T, Ko T, Sumida T, Toko H, Naito AT, Takeda N, Tobita T et al: High-throughput single-molecule RNA imaging analysis reveals heterogeneous responses of cardiomyocytes to hemodynamic overload. J Mol Cell Cardiol 2019, 128:77-89.
Friedman CE, Nguyen Q, Lukowski SW, Helfer A, Chiu HS, Miklas J, Levy S, Suo S, Han JJ, Osteil P et al: Single-Cell Transcriptomic Analysis of Cardiac Differentiation from Human PSCs Reveals HOPX-Dependent Cardiomyocyte Maturation. Cell Stem Cell 2018, 23(4):586-598.e588.
Yekelchyk M, Guenther S, Preussner J, Braun T: Mono- and multi-nucleated ventricular cardiomyocytes constitute a transcriptionally homogenous cell population. Basic Res Cardiol 2019, 114(5):36.
Bhullar RP, Clough RR, Kanungo J, Elsaraj SM, Grujic O: Ral-GTPase interacts with the beta1 subunit of Na+/K+-ATPase and is activated upon inhibition of the Na+/K+ pump. Can J Physiol Pharmacol 2007, 85(3-4):444-454.
Barwe SP, Jordan MC, Skay A, Inge L, Rajasekaran SA, Wolle D, Johnson CL, Neco P, Fang K, Rozengurt N et al: Dysfunction of ouabain-induced cardiac contractility in mice with heart-specific ablation of Na,K-ATPase beta1-subunit. J Mol Cell Cardiol 2009, 47(4):552-560.
Bass-Stringer S, Ooi JYY, McMullen JR: Clusterin is regulated by IGF1-PI3K signaling in the heart: implications for biomarker and drug target discovery, and cardiotoxicity. Arch Toxicol 2020.
Ullrich M, Aßmus B, Augustin AM, Häbich H, Abeßer M, Martin Machado J, Werner F, Erkens R, Arias-Loza AP, Umbenhauer S et al: SPRED2 deficiency elicits cardiac arrhythmias and premature death via impaired autophagy. J Mol Cell Cardiol 2019, 129:13-26.
Kakimoto Y, Sasaki A, Niioka M, Kawabe N, Osawa M: Myocardial cathepsin D is downregulated in sudden cardiac death. PLoS One 2020, 15(3):e0230375.
Fujita M, Mitsuhashi H, Isogai S, Nakata T, Kawakami A, Nonaka I, Noguchi S, Hayashi YK, Nishino I, Kudo A: Filamin C plays an essential role in the maintenance of the structural integrity of cardiac and skeletal muscles, revealed by the medaka mutant zacro. Dev Biol 2012, 361(1):79-89.
Kebir S, Orfanos Z, Schuld J, Linhart M, Lamberz C, van der Ven PF, Schrickel J, Kirfel G, Fürst DO, Meyer R: Sarcomeric lesions and remodeling proximal to intercalated disks in overload-induced cardiac hypertrophy. Exp Cell Res 2016, 348(1):95-105.
Shinde AV, Su Y, Palanski BA, Fujikura K, Garcia MJ, Frangogiannis NG: Pharmacologic inhibition of the enzymatic effects of tissue transglutaminase reduces cardiac fibrosis and attenuates cardiomyocyte hypertrophy following pressure overload. J Mol Cell Cardiol 2018, 117:36-48.
Oh YJ, Pau VC, Steppan J, Sikka G, Bead VR, Nyhan D, Levine BD, Berkowitz DE, Santhanam L: Role of tissue transglutaminase in age-associated ventricular stiffness. Amino Acids 2017, 49(3):695-704.
You J, Wu J, Zhang Q, Ye Y, Wang S, Huang J, Liu H, Wang X, Zhang W, Bu L et al: Differential cardiac hypertrophy and signaling pathways in pressure versus volume overload. Am J Physiol Heart Circ Physiol 2018, 314(3):H552-h562.
Fan C, Li Y, Yang H, Cui Y, Wang H, Zhou H, Zhang J, Du B, Zhai Q, Wu D et al: Tamarixetin protects against cardiac hypertrophy via inhibiting NFAT and AKT pathway. J Mol Histol 2019, 50(4):343-354.
Gesmundo I, Miragoli M, Carullo P, Trovato L, Larcher V, Di Pasquale E, Brancaccio M, Mazzola M, Villanova T, Sorge M et al: Growth hormone-releasing hormone attenuates cardiac hypertrophy and improves heart function in pressure overload-induced heart failure. Proc Natl Acad Sci U S A 2017, 114(45):12033-12038.
Anderson ME, Brown JH, Bers DM: CaMKII in myocardial hypertrophy and heart failure. J Mol Cell Cardiol 2011, 51(4):468-473.
Gibb AA, Hill BG: Metabolic Coordination of Physiological and Pathological Cardiac Remodeling. Circ Res 2018, 123(1):107-128.
Chu W, Wan L, Zhao D, Qu X, Cai F, Huo R, Wang N, Zhu J, Zhang C, Zheng F et al: Mild hypoxia-induced cardiomyocyte hypertrophy via up-regulation of HIF-1α-mediated TRPC signalling. J Cell Mol Med 2012, 16(9):2022-2034.
Zhu ZY, Gao T, Huang Y, Xue J, Xie ML: Apigenin ameliorates hypertension-induced cardiac hypertrophy and down-regulates cardiac hypoxia inducible factor-lα in rats. Food Funct 2016, 7(4):1992-1998.

TableS5.pdf
Table S5. Hub gene detection for TAC development.
TableS4.pdf
Table S4. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment terms for TAC development.
TableS3.pdf
Table S3. The Gene Ontology (GO) enrichment terms for TAC development.
TableS2.pdf
Table S2. Hub markers detection for the different cell cltusers in cardiac hypertrophy development.
TableS1.pdf
Table S1. The GSEA enrichment terms for transverse aortic constriction (TAC) development.

Download PDF

Version 1

posted

You are reading this latest preprint version

Single-cell RNA-sequencing Unravels Heterogeneity of Cardiomyocytes and Signaling Pathways in Pressure Overload

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

Background

Materials and methods

Results

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1