Single cardiomyocyte nuclear transcriptomes reveal a lincRNA-regulated de-differentiation and cell cycle stress-response in vivo

Abstract

Cardiac regeneration may revolutionize treatment for heart failure but endogenous progenitor-derived cardiomyocytes in the adult mammalian heart are few and pre-existing adult cardiomyocytes divide only at very low rates. Although candidate genes that control cardiomyocyte cell cycle re-entry have been implicated, expression heterogeneity in the cardiomyocyte stress-response has never been explored. Here, we show by single nuclear RNA-sequencing of cardiomyocytes from both mouse and human failing, and non-failing adult hearts that sub-populations of cardiomyocytes upregulate cell cycle activators and inhibitors consequent to the stress-response in vivo. We characterize these subgroups by weighted gene co-expression network analysis and discover long intergenic non-coding RNAs (lincRNA) as key nodal regulators. KD of nodal lincRNAs affects expression levels of genes related to dedifferentiation and cell cycle, within the same gene regulatory network. Our study reveals that sub-populations of adult cardiomyocytes may have a unique endogenous potential for cardiac regeneration in vivo.Adult mammalian cardiomyocytes are predominantly binucleated and unable to divide. Using single nuclear RNA-sequencing of cardiomyocytes from mouse and human failing and non-failing adult hearts, See et al. show that some cardiomyocytes respond to stress by dedifferentiation and cell cycle re-entry regulated by lncRNAs.

Fig. 1

Single nuclear RNA-seq reveals CM heterogeneity. a, b Core cardiac genes that are most highly expressed in every CM nucleus (a) exhibit high expression with low coefficient of variation (b). c Highly expressed genes in TAC nuclei have higher penetrance than highly expressed genes in Sham nuclei. Spearman’s rank correlation (r = 0.90, p < 2.2e⁻¹⁶) shows good correlation between average expression level and penetrance. d Density distribution of correlation shows higher correlation in TAC nuclei than in Sham nuclei. p < 2.2e⁻¹⁶ from Mann–Whitney U test. e, f Unsupervised hierarchical clustering e and PCA f of single nuclear RNA-seq of CM reveal that CM nuclei accurately segregate into clusters specific to Sham or TAC subgroups (subgroup a, b) and is replicated across biological repeats (Rep) f. g Ranked Spearman correlation plot shows higher correlation in TAC nuclei than in Sham nuclei, which is replicated across biological repeats (Rep)

Fig. 2

LincRNAs in nodal hubs of gene regulatory networks. a, b WGCNA identifies three distinct gene modules (Healthy, Disease 1 and Disease 2) (a) in Sham and TAC nuclei that represent expression signatures of specific Sham or TAC nuclear subgroups (b). c–e WGCNA reveals candidate lincRNAs in nodal hubs bearing the highest connectivity with other genes within the gene regulatory network modules. Gas5 and Sghrt are in nodal hubs within disease module 2 (e) and highly correlated with expression of other genes in the network such as Nppa, Dstn, Ccng1, and Ccnd2. Size and color of bubbles represent strength and significance of connectivity. Key enriched gene ontology (GO) terms are listed for each module (p < 0.05 Fisher’s exact test). f–h Scatterplots showing the expression of genes from the 3 gene modules at the single-nuclear level (f), at pooled nuclei level (g) and matched bulk left ventricle tissue RNA-seq (h). i Significant differential expression of genes from the three gene modules between Sham and TAC samples is detected only by single nuclear RNA-seq, and not by pooled nuclei or bulk tissue RNA-seq

Fig. 3

Human single cardiomyocyte nuclear RNA-seq. a Core cardiac genes in human CMs are similar to mouse. b–d Unsupervised hierarchical clustering (b), PCA (c) and Spearman correlation analysis (d) produced 2 distinct subgroups in each of control and dilated cardiomyopathy (DCM) nuclei. e Density distribution of correlation shows narrower distribution for DCM nuclei compared to control. P value from Mann–Whitney U test. f WGCNA identifies gene modules (healthy 1, healthy 2, disease 1, and disease 2) that are specific for DCM or control nuclear subgroups. g, h Classifiers from human gene modules show differential expression at single nuclear level (g), but not in matched bulk left ventricle RNA-seq (h)

Fig. 4

Quadrant analyses reveal sub-populations of CM. a–c Quadrant analysis for Proliferation vs. Negative regulators of proliferation genes identifies increased co-expression in individual TAC nuclei (Q2; 44.4%; p = 3.237e⁻⁰⁷ Fisher’s exact test), only detectable by single nuclear RNA-seq (a), and not in pooled nuclei (b) or matched bulk left ventricle RNA-seq (c). Inset: histogram of nuclei distributed across quadrants. Blue represents Sham and red represents TAC nuclei. d–f Quadrant analysis for cardiac progenitor vs. cardiac transcription factor gene expression shows increased co-expression upon TAC stress in single CM nuclei d (Q2; 42.9%; p = 2.548e⁻⁰⁵ Fisher’s exact test), again not detectable in pooled nuclei or bulk tissue RNA-seq (e, f). g–i Increased co-expression of fetal reprogramming genes and dedifferentiation markers under TAC stress only detected in single nuclear RNA-seq (g) (Q2; 58.73%; p = 0.001371 Fisher’s exact test) and not in non-single approaches (h, i). j High co-expression of cardiac progenitors, cardiac transcription factors, dedifferentiation, proliferation, and negative proliferation markers is confined to single nuclear TAC samples in Q2 and Q4. k, l Single molecule RNA FISH shows Sca1 upregulation and co-expression of Tnnt2 in isolated adult mouse CMs from TAC hearts l compared to Sham k. Number of Sca1 ⁺ Sham CMs: 5/13; Sca1 ⁺ TAC CMs: 38/55; all together from 2 Sham and 3 TAC biological replicates. m, n Immunofluorescence confirms increase in cell-surface SCA1 protein expression in TAC CMs (n) compared to Sham (m). Number of SCA1⁺ Sham CMs: 8/23; SCA1⁺ TAC CMs: 43/66; all together from 2 Sham and 3 TAC biological replicates. Scale bar represents 20 μm

Fig. 5

LINCM expression validated by single molecule RNA FISH. a Single nuclear RNA-seq identifies 141 novel lincRNAs in nuclei of CMs (LINCMs) that are not in current public databases (ENSEMBL and NONCODE) nor published cardiac transcriptome data sets. b Single nuclear RNA-seq identifies LINCMs that are not detectable in matched left ventricle bulk tissue RNA-seq, explained by the dilution of reads in cytoplasmic mRNA pool. c Active H3K27Ac enhancer chromatin regions proximal to LINCMs are enriched in MEF2 transcription factor binding motif and functionally annotated by GREAT analysis to have cardiac expression and phenotypes. d Sites of active transcription demonstrated by co-localization of exonic and intronic probes (asterisk) in nucleus. Scale bar represents 5 μm. e–m Single molecule RNA FISH validates the expression of LINCMs in isolated adult mouse CMs. n–q Positive controls for highly abundant core genes Tpm1, Tnnt2, Myl2, and Malat1. r, s Negative controls with no-probe control (NPC) (r) and sense probe (s) to confirm signal specificity. Scale bar represents 10 μm. t, u Gas5 is upregulated in TAC CM and co-localizes with perinuclear Nppa transcripts. v, w Sghrt is upregulated and localizes to the cytoplasm of TAC CM. x, y LINCM5 is downregulated in TAC CM. Scale bar represents 10 μm

Fig. 6

Gas5 and Sghrt regulate co-regulated network gene expression. a Strategy to KD Gas5 or Sghrt independently in isolated adult TAC CMs in vitro. Cartoon from openclipart.org, under a CC0 1.0 Universal license. b, c In vitro KD of Gas5 or Sghrt in adult mouse CMs using GapmeRs is efficient and reproducible across biological replicates. N = 5 biological replicates. d–g Gas5 KD in TAC CMs results in significant reduction of Nppa, Dstn, Ccng1, and Ccnd2 expression. Sghrt KD in TAC CMs results in significant increase in Ccng1 and reduction in Ccnd2 expression. N = 5 biological replicates. Data represented as mean ± s.e.m. *p < 0.05, **p < 0.01, and ***p < 0.001