Single-nucleus transcriptomic survey of cell diversity and functional maturation in postnatal mammalian hearts.

A fundamental challenge in understanding cardiac biology and disease is that the remarkable heterogeneity in cell type composition and functional states have not been well characterized at single-cell resolution in maturing and diseased mammalian hearts. Massively parallel single-nucleus RNA sequencing (snRNA-seq) has emerged as a powerful tool to address these questions by interrogating the transcriptome of tens of thousands of nuclei isolated from fresh or frozen tissues. snRNA-seq overcomes the technical challenge of isolating intact single cells from complex tissues, including the maturing mammalian hearts; reduces biased recovery of easily dissociated cell types; and minimizes aberrant gene expression during the whole-cell dissociation. Here we applied sNucDrop-seq, a droplet microfluidics-based massively parallel snRNA-seq method, to investigate the transcriptional landscape of postnatal maturing mouse hearts in both healthy and disease states. By profiling the transcriptome of nearly 20,000 nuclei, we identified major and rare cardiac cell types and revealed significant heterogeneity of cardiomyocytes, fibroblasts, and endothelial cells in postnatal developing hearts. When applied to a mouse model of pediatric mitochondrial cardiomyopathy, we uncovered profound cell type-specific modifications of the cardiac transcriptional landscape at single-nucleus resolution, including changes of subtype composition, maturation states, and functional remodeling of each cell type. Furthermore, we employed sNucDrop-seq to decipher the cardiac cell type-specific gene regulatory network (GRN) of GDF15, a heart-derived hormone and clinically important diagnostic biomarker of heart disease. Together, our results present a rich resource for studying cardiac biology and provide new insights into heart disease using an approach broadly applicable to many fields of biomedicine.

Significant genetic and functional heterogeneity exist among individual cells and contribute to overall organ development, physiology, and plasticity (Altschuler and Wu 2010). For example, there are many different cell types, including cardiomyocytes, fibroblasts, endothelial cells (ECs), epicardial cells, blood cells, and others, in both developing and adult hearts, with significant anatomical, genetic, and functional heterogeneity among the same cell type (Paige et al. 2015; Potente and Makinen 2017; Tallquist and Molkentin 2017). Of particular interest are early postnatal hearts, in which the cell type composition and metabolic states undergo significant changes and remodeling for functional maturation. However, the postnatal period is significantly underexplored compared with the embryonic period or adulthood. In addition, cell type-specific metabolic and transcriptional remodeling often occurs in the disease state. Many forms of heart disease, such as mitochondrial cardiomyopathy, a disease of heart muscle due to primary mitochondrial dysfunction, are well recognized to exhibit significant heterogeneity in disease genetics, symptoms, and outcomes (Meyers et al. 2013; Burke et al. 2016b; Ritterhoff and Tian 2017). However, the precise cell type composition and molecular underpinnings of such cellular and functional heterogeneity in maturing and diseased postnatal hearts remain poorly understood.

Single-cell RNA sequencing (RNA-seq) approaches empower genome-wide gene expression analysis at the single-cell resolution and have provided novel insights into many important biological questions, such as cellular identity and heterogeneity (Tanay and Regev 2017; Woodworth et al. 2017). Recent technological advances in droplet microfluidics have enabled massively parallel single-cell transcriptome analysis of tens of thousands of cells (Klein et al. 2015; Macosko et al. 2015) and nuclei (Lake et al. 2016; Habib et al. 2017; Hu et al. 2017). Although single-nucleus RNA-seq (snRNA-seq) often detects a lower number of RNA than single-cell RNA-seq, as it excludes transcripts outside the nucleus, it holds several technical advantages when applied to tissues that are difficult to dissociate, such as mammalian hearts. First, snRNA-seq bypasses preparation of intact single-cell suspension and thus can be used to study frozen and archived primary tissues. Second, snRNA-seq can minimize the bias of recovering only cells that are easily dissociated or cells of a certain size optimal for platform-specific cell capture and reduce the impact of aberrant transcriptional changes induced by enzymatic dissociation. To date, massively parallel snRNA-seq has been used only to study brain cell type compositions and functional states (Lake et al. 2016; Habib et al. 2017; Hu et al. 2017). Further applications of massively parallel snRNA-seq to study cardiac biology and disease will provide novel insights into physiologically and clinically significant questions that cannot be satisfactorily addressed with population-based genome-wide profiling methods.

Here we applied our recently developed massively parallel snRNA-seq method, sNucDrop-seq (Hu et al. 2017), to investigate the transcriptional landscape of postnatal mouse hearts in both healthy and disease states of mitochondrial cardiomyopathy. Our transcriptome analysis of nearly 20,000 nuclei offered new insights into several key questions in cardiac biology: cell type composition, and functional heterogeneity in normal postnatal hearts and changes of transcriptional landscape of each cell type in disease state. We also used our single-nucleus transcriptome data set to reveal cell type-specific gene regulatory networks (GRNs) exemplified by heart disease-associated Gdf15 transcription. Our approach is applicable to study similar questions in many areas of biology and disease.

Results

sNucDrop-seq for single-nucleus transcriptome analysis of postnatal mouse hearts

We optimized a mouse heart nucleus isolation protocol based on sucrose gradient ultracentrifugation that helps minimize cytoplasmic contamination and protect nucleus integrity (Supplemental Fig. S1A; Hu et al. 2017). We performed sNucDrop-seq in normal developing postnatal hearts as well as hearts from a mouse model of pediatric mitochondrial cardiomyopathy. In this model, cardiac genetic inactivation of two transcription factors essential for normal cardiac metabolism and function (estrogen-related receptor α [ERRα] and ERRγ) results in rapid postnatal development of dilated mitochondrial cardiomyopathy, heart failure, and death within a month of birth (Wang et al. 2015). ERRα and ERRγ directly regulate expression of hundreds of genes important in mitochondrial fatty acid oxidation and oxidative phosphorylation (OxPhos) as well as cardiac contraction and conduction (Alaynick et al. 2007; Dufour et al. 2007; Huss et al. 2007; Wang et al. 2015). Cardiac ERRα/γ knockout (referred to here as knockout) mouse hearts exhibited loss of mitochondrial structure and function as well as defects of myocardial contraction and conduction, accompanied by significantly reduced expression of mitochondrial and cardiac function genes (Wang et al. 2015). To optimize and validate the sNucDrop-seq assay for postnatal heart tissues, we performed sNucDrop-seq analysis of dissected ventricles from control and knockout mice (n = 3 littermate pairs) of 9–10 d of age—an early stage of disease development in knockout, when significant gene expression and functional changes could be readily detected (Wang et al. 2015, 2017).

We performed sNucDrop-seq of both freshly isolated (control 1 and knockout 1) and frozen (control 2 and 3 and knockout 2 and 3) heart samples and obtained highly concordant results within the same genotype (Supplemental Fig. S1B,C). Overall, 78% of reads aligned to genomes, among which 77% mapped to exons, 16% mapped to introns, and 7% mapped to intergenic regions. This relatively lower percentage of reads mapped to the intronic region in the nuclear transcriptomic profiles of heart samples (compared with ∼50% intronic reads in mouse brains) (Hu et al. 2017) suggests that the relative composition of nascent transcripts varies significantly among cell types and organs. After quality filtering (>500 genes detected per nucleus), >15,000 nuclei were retained from three pairs of control and knockout littermates (Supplemental Table S1) for further analysis (7760 nuclei for control and 7323 nuclei for knockout). We obtained similar numbers and distributions of transcripts and genes per nucleus between samples (Supplemental Fig. S1B; Supplemental Table S1). In addition, sNucDrop-seq results showed high concordance when compared with bulk RNA-seq from control and knockout hearts (Supplemental Fig. S1D), further validating the sNucDrop-seq approach. sNucDrop-seq also provided additional, previously inaccessible insights into these transcriptional changes at single-nucleus resolution: Differential gene expression changes (e.g., Nppa and Cox7a1) between control and knockout hearts could result from changes in both the percentage of positive-expressing cells and average gene expression per nucleus (Supplemental Fig. S1E).

Unbiased cell type identification and study of cellular heterogeneity in healthy postnatal hearts

We assigned 7760 nuclei from healthy hearts into 14 distinct clusters using principal component analysis (PCA) dimension reduction followed by graph-based clustering, which was visualized by spectral t-distributed stochastic neighbor embedding (tSNE) (Fig. 1A). Numbers of detected transcripts (unique molecular identifiers [UMIs]) and genes were comparable among all clusters (Supplemental Fig. S2A,B). Each cluster contained nuclei from all heart samples, indicating overall reproducible transcriptional identities across biological replicates. Each of these 14 clusters displayed distinct nuclear transcriptional signatures, including both protein-coding and noncoding RNAs (Fig. 1B; Supplemental Fig. S2C,D; Supplemental Table S2). Based on pathway analysis and cell type-specific markers (Fig. 1C,D), we identified all of the major cardiac cell types—including cardiomyocytes, fibroblasts, and ECs—that express cell type-specific marker genes and together constituted >90% of all nuclei. Cardiomyocytes could be divided into two major groups based on mitochondrial content, muscle fiber type, and key markers: less mature or developing cardiomyocytes (dCMs) with fewer mitochondria and positive for immature cardiomyoctye markers such as Myocd (also known as myocardin) and more mature cardiomyocytes (mCMs) with abundant mitochondria and positive for muscle fiber markers such as Actc1 (also known as cardiac α-actin). Importantly, the relative cell type composition uncovered by sNucDrop-seq agreed well with the results defined by orthogonal approaches, including immunohistochemistry, FACS, and lineage tracing (Banerjee et al. 2007; Doppler et al. 2017). For instance, it was reported previously that 15-d-old (postnatal day 15 [P15]) mouse hearts contained 63% cardiomyocytes and 18% fibroblasts (Banerjee et al. 2007); we identified 59% cardiomyocytes and 19% fibroblasts in P10 mouse hearts.

Figure 1.

Unbiased cell type identification in the postnatal heart. (A) tSNE plot of 14 clusters of a P10 control heart. Cell identity and percentage are labeled. (B) Heat map showing clustering and gene expression of the 14 clusters. (C,D) Violin plot (C) and feature plot (D) illustrating the expression patterns of selected marker genes of each cluster. (pCM) Proliferating cardiomyocytes; (LEC) lymphatic ECs; (Fb) fibroblasts; (Epi) epicardial cells; (BC) blood cells; (PC/SMC) pericytes/smooth muscle cells.

Importantly, sNucDrop-seq revealed significant cellular heterogeneity within each major cell type. We were able to identify three populations of dCMs as well as two populations of mCMs, fibroblasts, and ECs. To determine whether these are distinct clusters, we performed differential expression tests to identify cell subtype-specific marker genes, which revealed subtle but significant differences between cell subtypes (Supplemental Table S3). Among dCM subpopulations, dCM3 was more mature, as it expressed more mitochondrial genes (e.g., Cox6a2, Ndufa5, and Atp5j), while dCM1 expressed fibroblast-enriched markers (e.g., Col1a2, Col3a1, and Dcn). Interestingly, some mCMs, such as those in the mCM1, also expressed these fibroblast-enriched markers. There was also a small population of presumably proliferating cardiomyocytes (pCMs) that expressed cell cycle genes (e.g., Mki67, Cenpp, and Kif15). Over 50% of these cells also expressed Gata4 and Myocd, transcription factors critical for cardiomyocyte development (Fig. 1C,D). Indeed, Gata4+ or Myocd+ nuclei were significantly enriched only in dCMs (dCM1–3; odds ratio = 7.96; P < 2.2 × 10−16 by Fisher's exact test) but not in mCMs or nonmyocyte cells. Overall, these results reveal significant heterogeneity among dCMs, mCMs, and fibroblasts, with many subtypes that exhibited intermediate molecular signatures.

Through sNucDrop-seq analysis, we also identified two distinct EC populations not reported in previous single-cell RNA-seq studies of mouse embryonic (DeLaughter et al. 2016; Li et al. 2016) or adult (Gladka et al. 2018; Skelly et al. 2018) hearts. Analysis of subtype-enriched genes uncovered that EC2 expressed relatively more mitochondrial and muscle fiber genes (e.g., Cox6a2, Mybpc3, and Mhrt) (Supplemental Table S3). EC1-enriched genes included Npr3 as well as Tie1 and Tie2 (receptors for angiopoietins critical for angiogenesis) (Supplemental Fig. S3A). As a proof of principle for using our cardiac cell atlas as a resource to identify and characterize these EC subtypes, we used immunohistochemistry to show that these EC subtypes represent anatomically distinct populations of cardiac ECs (Supplemental Fig. S3B). While general EC marker PECAM1 stained all ECs in the endocardium, arteries, veins, and capillaries, EC1-specific NPR3 marked only some capillary and endocardium but few arteries or veins. We also used confocal microscopy to examine these ECs with double immunostaining. PECAM1 and NPR3 double staining revealed that many NPR3 cells are PECAM1-positive and therefore ECs (Supplemental Fig. S3C). PECAM1 and cardiac troponin I (TNNI3) double staining revealed that there are some cells that coexpress PECAM1 and TNNI3 (Supplemental Fig. S3D), suggesting that a subpopulation of ECs can express muscle fiber genes.

The sNucDrop-seq analysis also identified relatively rare cell types. These included pericytes or smooth muscle cells, blood cells, epicardial cells, and pCMs, which express related specific marker genes and together constitute 8.8% of the total nuclei (Fig. 1A–D). Historically, these rare cardiac cell types were often defined by a few cell markers. Our unbiased approach not only validated known markers (Fig. 1C,D; Supplemental Table S2) but also revealed transcriptomic signatures based on the expression patterns of hundreds of genes for these cell types. Based on our sNucDrop-seq data, we identified unbiased lists of highly cell type-specific genes that can be used to classify cardiac pericytes/smooth muscle cells, epicardial cells, and lymphatic ECs (LECs), respectively (Supplemental Table S4). For example, although Wt1 and Tbx18 are well-established epicardium markers, our study uncovered >30 additional high-confidence markers for epicardial cells, including both known (e.g., Upk3b and Upk1b) and novel (e.g., Il18r1 and Msln) genes (Supplemental Fig. S3E; Bochmann et al. 2010). Moreover, these new marker genes may reveal previously unrecognized biological functions of these less studied cell types; for example, mechanosensitive ion channel Piezo2 and genes important in neurodevelopment (e.g., Shank3 and Reln) were highly expressed and specifically enriched in LECs (Supplemental Fig. S3F). Indeed, 62% of the Flt4+ (Vegfc receptor, which is a marker for lymphatic vessel) cells coexpressed Piezo2, suggesting potential functions of some LECs in sensing blood pressure or other mechanical signals. This significantly expanded list of potential marker genes highlights the value of our study as a resource for future research investigations of these relatively rare cardiac cells.

Transcriptomic dynamics of cardiomyocyte maturation in postnatal hearts

To investigate transcriptome dynamics during postnatal heart development, we also performed sNucDrop-seq in P6 control mouse ventricles (n = 3 mice) and compared the results with those of P10 control mice. Using the same analysis settings, we obtained 4560 nuclei with similar numbers of transcripts and genes per nucleus among the three biological replicates (Supplemental Table S1). Ten distinct clusters were identified based on their transcriptomic signatures (Fig. 2A; Supplemental Table S5). All of the cell types identified in P10 hearts were found in P6 hearts except for LEC, likely because it is relatively rare (0.2% or 19 nuclei in P10) and therefore harder to detect, with fewer total nuclei studied (4560 nuclei in P6 vs. 7760 nuclei in P10). While two populations each of ECs and fibroblasts were also detected in P6 hearts, only one population each of dCMs and mCMs was found in P6 hearts. Pairwise comparison of all clusters between P6 and P10 hearts indicated that the cell types and subtypes identified in P6 hearts were counterparts of those in P10 hearts (Fig. 2B). Comparative analysis of cell type composition between P6 and P10 revealed that the biggest difference was the substantial decrease of pCMs and the corresponding increase of mCMs in P10 hearts, while the percentages of other cell types changed little (Fig. 2C). Notably, pCMs were much more abundant in P6 hearts than in P10 hearts (5.5% in P6 vs. 1.8% in P10). P6 and P10 pCM marker genes are highly overlapping, with 78% (101 out of 130) of P6 pCM marker genes also enriched in P10 pCMs. P6 and P10 pCMs exhibit similar gene signatures, dominated by essential cellular processes and pathways required for cell proliferation, such as mitotic cell cycle and chromosome segregation, and, in particular, cytokinesis for cell division (Fig. 2D). Ki67 (a pCM and proliferating cell marker) immunostaining confirmed the existence of these pCMs in mouse hearts, and their percentage declined during postnatal development (Fig. 2E), consistent with our sNucDrop-seq data. These results suggest that cardiomyocyte differentiation and/or maturation represents a major biological event during this period of postnatal heart development in mice (between P6 and P10).

Figure 2.

Transcriptomic dynamics of cardiomyocyte maturation in postnatal hearts. (A) tSNE plot of 10 clusters of a P6 control heart. Cell identities are labeled. (B) Pairwise comparison of all clusters between P6 and P10 control hearts. (C) Cell type compositions in P6 and P10 control hearts. (D) The top 10 enriched pathways in P6 and P10 pCMs. (E) Confocal microscopy shows that Ki67 staining (red) colocalizes with nucleus Hoechst staining (blue). (F, left panel) Workflow for PCA analysis. (Right panel) Plot of PC1 cell-loading scores (Y-axis). On the X-axis, cells are ordered by library size (largest to smallest) within each cell type (color-coded). (G) Heat map of the top loading genes. Nuclei are ordered by PC1. (Fb) Fibroblasts; (Epi) epicardial cells; (BC) blood cells; (PC/SMC) pericytes/smooth muscle cells.

To further dissect the developmental/transitional transcriptional changes within pCMs and dCMs, we performed PCA analysis on 5429 cardiomyocyte nuclei combined from |P6 (2205 nuclei) and P10 (3224 nuclei) using 1276 highly variable genes (Fig. 2F, left panel). The top PCA component, which accounted for 3.99% of total variation, separated all nuclei into overlapping stepwise clusters ordered by age (from P6 to P10) and cell type (from pCMs to dCMs), indicating the temporally progressive maturation of cardiomyocytes (Fig. 2F, right panel). We identified the top genes associated with PC1. These included cell cycle markers (i.e., Ckap2, Top2a, and Mki67) highly expressed in pCMs and muscle contraction genes (i.e., Ttn and Ryr2) highly enriched in dCMs (Fig. 2G). These results suggest that these pCM and dCM clusters may represent different stages along the postnatal cardiomyocyte maturation process.

Cell type-specific transcriptional remodeling of pediatric mitochondrial cardiomyopathy

We next applied sNucDrop-seq to investigate the cell type-specific transcriptional landscape in the cardiac ERRα/γ knockout mouse model of pediatric mitochondrial cardiomyopathy. Analysis of the 7323 nuclei from three P10 knockout hearts revealed 13 clusters (Fig. 3A) with largely comparable numbers of detected transcripts and genes for each cluster (Supplemental Fig. S4A). As in control hearts, each cluster contained a unique transcriptome, including protein-coding and noncoding RNA signatures (Fig. 3B,C; Supplemental Fig. S4B,C; Supplemental Table S6). Pathway analysis and marker genes revealed the identities of these clusters (Fig. 3A–C; Supplemental Table S6). Pairwise comparison of all control and knockout clusters showed that there was generally a one-to-one correlation of the same subtype of cardiac cells between control and knockout (Fig. 3D).

Figure 3.

Cell type identification and correlation in a mouse model of pediatric mitochondrial cardiomyopathy. (A) tSNE plot of 13 clusters of P10 knockout hearts. (B,C) Violin plot (B) and feature plot (C) illustrating the expression patterns of selected marker genes of each cluster. (D) Pairwise comparison of all clusters between P10 control and knockout hearts. (aFb) Activated fibroblasts; (Fb) fibroblasts; (Epi) epicardial cells; (BC) blood cells; (PC/SMC) pericytes/smooth muscle cells.

While knockout hearts retained most of the major and rare cell types except pCMs and LECs, their relative abundance exhibited cell type-specific changes (Fig. 4A). There remained significant heterogeneity among cardiomyocytes, fibroblasts, and ECs. However, each knockout cell type or subtype exhibited significant transcriptome-level changes (Fig. 4B). Subtypes of cardiomyocytes, fibroblasts, and ECs exhibited larger transcriptional changes than epicardial cells, pericytes/smooth muscle cells, or blood cells, suggesting that distinct cardiac cell types undergo differential levels of transcriptional remodeling during heart disease progression. Pathway enrichment analysis revealed that mitochondrial OxPhos and TCA cycle (e.g., Cox7a1, Ndufa1, Sdha), fatty acid degradation/oxidation (e.g., Acadm and Hadh), extracellular matrix associated with fibrosis (e.g., Post and Fbn1), and ribosome genes (e.g., Rplp1 and Rps3) were among the most affected pathways in knockout hearts (Fig. 4C,D). However, they also displayed cell type-specific and subtype-specific pathway changes. This suggests that individual cell types or subtypes adopted differential metabolic and functional adaptations to disease state.

Figure 4.

Cell type-specific transcriptional remodeling of pediatric mitochondrial cardiomyopathy. (A) Cell type compositions in P10 control and knockout hearts. (B) Volcano plot showing differentially expressed genes in correlated clusters of control and knockout hearts. (C) Cellular pathways significantly changed in knockout versus control hearts between correlating clusters. (D) Violin plot showing representative gene expression changes across different cell types. (E) Sirius Red stain of fibrosis in control and knockout hearts and quantification. n = 2 mice per genotype. Bar, 100 µm. (*) P < 0.05 by t-test. (Fb) Fibroblasts; (aFb) activated fibroblasts; (Epi) epicardial cells; (BC) blood cells; (PC/SMC) pericytes/smooth muscle cells; (BCAA) branched chain amino acid.

It has often been observed that mitochondrial functions and gene expression are decreased in mitochondrial cardiomyopathy and many other forms of heart disease, including our knockout mouse model (Huss and Kelly 2005; Abel and Doenst 2011; Harvey and Leinwand 2011; Wang et al. 2015). However, previous RNA-seq studies were based on whole hearts/chambers or focused only on cardiomyocytes that possess the most abundant mitochondria. Our sNucDrop-seq analysis revealed cell type-specific mitochondrial gene expression changes in diseased hearts for all major cell types at the single-nucleus resolution. Consistent with bulk RNA-seq and our previous report (Wang et al. 2015), most subtypes of dCMs and mCMs exhibited significantly decreased expression of mitochondrial OxPhos genes (Fig. 4C,D), reflecting the direct effects of loss of ERRα/γ in these cells. Reduced OxPhos gene expression was also observed in fibroblasts and the EC1 subtype of ECs in knockout hearts but not in blood cells, pericytes/smooth muscle cells, or epicardial cells. In contrast, expression of mitochondrial OxPhos genes was significantly increased in a putative population of activated fibroblasts, consistent with a metabolic profile of increased mitochondrial respiration reported recently in these cells (Bernard et al. 2015; Negmadjanov et al. 2015). These mitochondrial gene expression changes in noncardiomyocytes are likely secondary due to the progression of mitochondrial cardiomyopathy.

It is generally appreciated that the heart switches from mainly catabolizing fat in the healthy state to predominantly using glucose in the heart failure stage (Huss and Kelly 2005). In addition, defective branched chain amino acid (BCAA) catabolism has been identified recently as a metabolic hallmark of heart failure and contributes to disease development (Sun et al. 2016; Li et al. 2017). Our sNucDrop-seq results revealed that only mCMs, not dCMs or other cell types, showed a significant decrease in gene expression related to fatty acid oxidation, BCAA catabolism, and associated pyruvate metabolism (Fig. 4C,D; Li et al. 2017). This was in contrast to a more universal decrease in mitochondrial TCA and OxPhos. Therefore, the general concept of mitochondrial dysfunction in heart disease occurs in a more intricate cell type-specific manner: Defective fatty acid oxidation and BCAA catabolism are mCM-specific, but a decrease in mitochondrial OxPhos takes place in a broader spectrum of cardiac cell types. These results suggest that targeting metabolic remodeling for heart disease treatment needs to consider these cell type-specific mitochondrial pathways.

Increased expression of fibrosis-associated extracellular matrix genes (e.g., Postn and Fbn1) was observed across most cell types in knockout hearts compared with wild-type hearts (Fig. 4C,D). This suggests that cardiac fibrosis may occur in a much broader scope than generally believed: It affects the whole heart and is not limited to fibroblasts. Notably, we also observed an increase of ribosome genes expression in many cell types of diseased hearts. Increased protein synthesis or ribosome RNA abundance was reported previously in animal models of dilated or hypertrophic cardiomyopathy (Moroz 1967; Morgan et al. 1985; Burke et al. 2016a). Our results suggest that such a potentially compensatory response in protein synthesis occurs in mitochondrial cardiomyopathy and involves all major cardiac cell types.

Increased fibrosis is a hallmark of many forms of heart disease, including mitochondrial cardiomyopathy. Although knockout hearts retained a comparable percentage of total fibroblasts, a putative population of activated fibroblasts was identified based on known fibrosis marker genes in addition to the inactive fibroblast population. Activated fibroblasts exhibited transcriptional profiles of enriched fibrosis-related extracellular matrix as well as endoplasmic reticulum and mitochondrial OxPhos pathways, consistent with known phenotypes of activated fibroblasts (Supplemental Table S5; Bernard et al. 2015; Negmadjanov et al. 2015; Tallquist and Molkentin 2017). Indeed, a significantly increased level of fibrosis was observed in knockout hearts at this age, as revealed by Sirius Red staining (Fig. 4E). Detailed examination of the activated fibroblast transcriptional signature revealed that they contained few genes involved in cell proliferation (Supplemental Table S6). Considering that knockout hearts contained similar percentages of total fibroblasts (fibroblasts and activated fibroblasts) (Fig. 4A) and lacked distinct populations of proliferating cells, this suggests that the increased fibrosis in knockout hearts results from fibroblast activation rather than proliferation.

While knockout hearts retained two EC populations (EC1 and EC2) and the total percentage of ECs remained unchanged between control and knockout, the relative composition of these EC subtypes likely changed in response to disease states. EC1 showed a significant expansion (8.6% compared with 4.8% in control) in knockout hearts, which was validated by immunostaining of the EC1-specific marker NPR3 (Supplemental Fig. S5A). Expression of Tie1 and Tie2 remained higher in EC1 than in EC2 (Supplemental Fig. S5B), suggesting that EC1 abundance might increase in knockout hearts due to increased angiogenesis.

GRNs associated with the transcriptional activation of Gdf15

GRN refers to a network of molecular regulators that functionally interact to modulate a gene expression program critical to specific cellular states and functions (Wyrick and Young 2002). Single-cell transcriptome analysis enables the inference of cell type-specific or developmental stage-specific GRNs (Moignard et al. 2015; Wagner et al. 2016). Based on sNucDrop-seq results, we next investigated the GRNs underlying the transcriptional activation of Gdf15 in heart disease. GDF15 (also known as MIC-1, NAG-1, PLAB, or PTGFB) is a protein of the transforming growth factor β (TGFβ) family with pleiotropic functions (Unsicker et al. 2013). GDF15 has recently gathered strong interest for its pharmaceutical potential in treating obesity and its implication in cardiovascular disease. GDF15 suppresses appetite through activating its receptor, GFRAL, in the brain stem (Johnen et al. 2007; Emmerson et al. 2017; Hsu et al. 2017; Mullican et al. 2017; Xiong et al. 2017; Yang et al. 2017) and is currently being investigated as a potential therapeutic agent/target for obesity, metabolic disease, anorexia, and cachexia. In addition, many clinical studies have found that plasma GDF15 is significantly elevated in various heart diseases and is associated with increased morbidity and mortality; it is therefore used as an independent biomarker for heart disease, especially heart failure (Baggen et al. 2017; Wollert et al. 2017). However, how GDF15 level is regulated at the molecular level in heart disease remains unclear (Wollert et al. 2017).

We demonstrated recently that cardiac Gdf15 transcription is markedly activated during the progression of heart failure in our cardiac ERRα/γ knockout mouse model, and, in the pediatric period, GDF15 acts as a heart-derived hormone to inhibit liver growth hormone signaling, thereby coordinating cardiac function and postnatal body growth (Wang et al. 2017). Our previous quantitative PCR (qPCR) analysis showed that average Gdf15 mRNA level increased ∼10-fold in P10 ERRα/γ knockout mouse hearts (Wang et al. 2017). GDF15 protein immunostaining and quantification further showed that, in contrast to very few GDF15+ cells in control mouse hearts, ∼1.5% of cells in knockout mouse hearts express detectable GDF15 protein at 10 d of age (Fig. 5A,B). However, the exact cell types expressing Gdf15 and how Gdf15 transcription is activated in these specific cell types remain incompletely understood.

Figure 5.

Gdf15 expression in control and knockout hearts. (A) Representative pictures of GDF15 immunostaining in P10 hearts. (B) Quantification of Gdf15+ cells based on immunostaining in A. n = 2 mice per genotype. (C) Feature plot showing Gdf15+ nucleus distribution in P10 hearts. (D) Quantification of Gdf15+ nuclei in control and knockout hearts based on sNucDrop-seq. (E) Distribution of Gdf15+ nuclei in each cell type. (**) P < 2.2 × 10−16 by χ2 test in B and D.

sNucDrop-seq allowed us to address these questions in a cell type-specific manner. Overall, there was a 10-fold increase in the number of Gdf15+ nuclei reaching ∼3% in P10 knockout hearts (Fig. 5C,D). The majority (72%) of Gdf15+ nuclei was cardiomyocytes (mCMs and dCMs) (Fig. 5E; Supplemental Table S7). This is consistent with our previous immunostaining result showing that GDF15 protein colocalizes with TNNI3 (Wang et al. 2017). Small numbers of Gdf15+ nuclei were also found in other cell types, notably fibroblasts (both fibroblasts and activated fibroblasts) and ECs. This raises the question of cell type specificity of GRNs that control Gdf15 expression.

To define cell type-specific GRNs that activate Gdf15 transcription in knockout hearts, we reasoned that the core network components (i.e., transcription factors, cofactors, and chromatin modifiers) should be coexpressed in Gdf15+ cells. We therefore searched for transcription factor, cofactor, and chromatin-modifying genes that (1) were statistically enriched in Gdf15+ nuclei compared with Gdf15− nuclei in either control or knockout hearts in the same cell type and (2) showed strong positive correlation of expression with Gdf15 in the same cell type. This analysis resulted in a list of candidate genes that constitute cell type-specific GRNs of Gdf15 (Fig. 6A; Supplemental Table S8). Expression of Gdf15 and these cell type-specific candidate GRN genes are highly correlated in cell types where both are present (Supplemental Fig. S6). For example, Gata4 was revealed as a candidate cardiomyocyte and fibroblast GRN component, and Gata4 and Gdf15 expression was positively correlated in both cardiomyocytes and fibroblasts. Notably, Gata4 is well known for its essential role in embryonic heart development, and mutations of the Gata4 gene are known to cause congenital heart disease in humans. In addition, the transcriptional activity of Gata4 is enhanced by pressure overload, angiotensin, and other signals that induce adult cardiomyopathy or heart failure through post-translational modification or increased expression (Oka et al. 2007). We next screened the mouse Gdf15 promoter (within 2 kb upstream of the transcription start site) and candidate enhancer regions (high H3K4me1 and H3K27ac histone marks close to the Gdf15 gene locus) and found that there were canonical binding motifs for almost all of the transcription factors of the candidate Gdf15 GRNs (Supplemental Table S9). Therefore, these transcription factors may directly bind to the Gdf15 promoter or enhancer to activate its transcription. Further analysis revealed the subset of genes with extensive known physical and functional interactions. Together, these genes form the nodes of the cell type-specific cardiac Gdf15 GRNs (Fig. 6B).

Figure 6.

Cardiac cell type-specific GRNs of Gdf15. (A) Venn diagram showing the number of transcriptional regulators that constitute cell type-specific GRNs. (B) Cardiac Gdf15 GRNs that include cell type-specific transcriptional regulators. Only genes of the key node status are shown. (C) Adenoviral-mediated overexpression of Gata4 induces Gdf15 expression in HL1 cells. (D) Gata4 activates Gdf15 expression in wild-type mouse hearts. (E) Increased binding of GATA4 to the mouse Gdf15 promoter by chromatin immunoprecipitation (ChIP). (F) Luciferase reporter assay showing that GATA4 activates the mouse Gdf15 promoter but not one with a mutant GATA4-binding site. The fibroblasts in A and B include both fibroblast and activated fibroblast cell populations. (*) P < 0.05; (**) P < 0.01 by t-test.

As a proof of principle, we next validated the biological significance of these findings by focusing on Gata4. We found that adenoviral-mediated Gata4 overexpression induced Gdf15 expression in HL-1 cardiomyocytes in vitro (Fig. 6C) and, importantly, in wild-type postnatal mouse hearts in vivo (Fig. 6D). Chromatin immunoprecipitation (ChIP) studies showed that GATA4 bound to a canonical GATA site in the mouse Gdf15 promoter (Fig. 6E). GATA4 increased activity of the Gdf15 promoter but not one with the GATA site mutated (Fig. 6F). Together, these results reveal the molecular mechanism of Gdf15 induction in heart disease and demonstrate the utility of snRNA-seq analysis in inferring cell type-specific GRNs.

Discussion

Here we applied massively parallel snRNA-seq to interrogate the transcriptional landscape in postnatal mouse hearts associated with either healthy or mitochondrial cardiomyopathy disease states. By resolving in vivo cellular heterogeneity in hearts with sNucDrop-seq, our study provided novel insights into several areas of cardiac biology that were underexplored. First, we showed that our unbiased approach could identify almost all cell types (abundant or rare), determine cell type composition of the heart, and resolve significant cellular heterogeneity associated with different developmental stages or functions in healthy postnatal hearts. Second, we investigated a pediatric mitochondrial cardiomyopathy mouse model, revealing distinct changes in organ composition and cell type-specific transcriptional remodeling with functional implications. Finally, we used the single-nucleus transcriptome results to uncover cell type-specific GRNs associated with transcriptional activation of Gdf15 in heart disease. Importantly, our approach is generally applicable to studying many areas of cardiac biology and disease.

The high throughput of sNucDrop-seq enabled transcriptomic analysis of ∼20,000 cardiac nuclei isolated from maturing postnatal hearts. Recent single-cell RNA-seq methods based on Fluidigm or flow cytometry platforms studied single heart cell transcriptome ranging from tens to ∼1000 cells per age/developmental stage during mouse embryonic heart development (DeLaughter et al. 2016; Li et al. 2016; Gladka et al. 2018). Another report also used the Fluidigm platform for snRNA-seq analysis of isolated adult cardiomyocytes with 13–64 nuclei per sample (See et al. 2017). Our study complemented these previous studies by focusing on the early postnatal period, when the heart is rapidly maturing. In addition, sNucDrop-seq overcame the technical challenges associated with isolating intact single cells from mammalian hearts, thus empowering many new insights: (1) We identified several relatively rare cell types, including blood cells, pericytes/smooth muscle cells, LECs, epicardial cells, and pCMs, and determined their relative abundance in postnatal hearts. (2) We discovered significant heterogeneity within cardiomyocytes, fibroblasts, and ECs. Previous single-cell RNA-seq studies identified in mouse embryonic hearts heterogeneous cardiomyocyte populations, with one cardiomyocyte subpopulation that expressed fibroblast genes (DeLaughter et al. 2016). Our results suggest an even broader degree of cellular heterogeneity, possibly due to higher throughput and the nature of the postnatal stage. (3) We revealed the transcriptome-wide changes of all cell types and subtypes in the disease state of mitochondrial cardiomyopathy, which exhibited significant subtype-specific remodeling.

The transcriptome-wide insights obtained by sNucDrop-seq also informed new insights into cardiac biology in both healthy and disease states. For example, we uncovered numerous new marker genes for rare cell types or heterogeneous cell subtypes. Compared with previous methods that depended on a few cell markers, our method is likely more accurate because it is based on transcriptome-level information. These new gene signatures are potentially associated with previously unrecognized functions of these cells and will facilitate future studies. We also resolved distinct transcriptional signatures between cell subtypes, such as EC1 and EC2, and validated the result by immunostaining. Our results presented here as a resource should assist future endeavors to investigate potential functional, anatomical, or developmental differences between cardiac cell subtypes. sNucDrop-seq also has the ability to capture ongoing transcriptional events compared with steady-state transcriptome analysis by whole-cell RNA-seq.

We identified populations of presumable pCMs in both P6 and P10 control hearts. It is thought that mouse hearts continues to grow by means of cardiomyocyte proliferation in the early postnatal period (Eschenhagen et al. 2017). These two pCM populations are enriched for cardiomyocyte development genes such as Gata4+ or Myocd+, suggesting that they are cardiomyocytes. In addition, these two pCM populations share the same transcriptional signature of cell proliferation: mitosis, G2/M transition, chromosome segregation, and, in particular, cytokinesis for cell division (Fig. 2D). Ki67 immunostaining corroborated the presence of these pCMs in mouse hearts (Fig. 2E). Furthermore, the abundance of these pCMs declines from P6 (5.5%) to P10 (1.8%), consistent with the notion that mouse cardiomyocytes progressively lose their ability to proliferate shortly after birth. However, it remains possible that at least some of these cells are simply undergoing DNA synthesis and mitosis without cell division, such as multinucleation. Therefore, caution is warranted regarding the exact nature of these cells, and additional studies are needed to definitively address the question of to what extent these pCMs are proliferating.

We were able to reveal cell type-specific and subtype-specific transcriptional remodeling in a mouse model of mitochondrial cardiomyopathy. Overall, cardiomyocytes, fibroblasts, and ECs exhibited profound metabolic and functional changes, while fewer changes were observed in other cardiac cell types. A common theme to most cells is the down-regulation of the mitochondrial OxPhos pathway (but more limited for fatty acid oxidation and BCAA catabolism pathways) and an increase of ribosome related genes. Increased fibrosis gene expression was also observed in many nonfibroblast cells, suggesting that tissue fibrosis might involve a broader range of cell types than simply fibroblast activation—at least at the transcriptional level. Intriguingly, we could not find signs of a distinct pCM population in knockout hearts any more. It is possible that mitochondrial cardiomyopathy had more profound impacts on cardiac health by impairing postnatal cardiomyocyte proliferation capacity (Naqvi et al. 2014; Eschenhagen et al. 2017).

The sNucDrop-seq data set allowed us to deduce candidate cell type-specific GRNs that activate Gdf15 transcription in heart disease, and we experimentally validated key nodes of these GRNs. Intriguingly, different cell types seem to use both overlapping (i.e., Gata4 for cardiomyocytes and fibroblasts) and cell type-specific (i.e., Gata6 for cardiomyocytes) transcription factors to regulate Gdf15 expression. One explanation is that these transcription factors exhibit cell type-specific expression patterns or their expression is induced only in specific cell types during heart disease progression. This raises the question of how these upstream regulators are transcriptionally activated in a cell type-specific manner (i.e., cell type-specific Gata4 GRNs). Based on our current study of Gdf15 GRNs, we anticipate that future time-course snRNA-seq studies covering multiple stages of heart disease progression will help address this question and shed light on the molecular underpinnings of heart disease progression. In summary, while future studies are needed to continue validating our current findings, our results demonstrate the value of large-scale snRNA-seq analysis in studying multiple areas of cardiac biology and disease.

Materials and methods

Animal studies

All animal studies were approved by and performed under the guidelines of the Institutional Animal Care and Use Committee of the Children's Hospital of Philadelphia (CHOP). All mice were in the C57BL6/J background. Cardiac ERRα/γ knockout mice were described previously (Wang et al. 2015).

snRNA-seq and data analysis

We performed sNucDrop-seq and data analysis as described previously (Hu et al. 2017) with optimizations for the heart tissue. Detailed methods and analysis are in the Supplemental Material. The data were deposited in the Gene Expression Omnibus (GEO) database under accession number GSE118545.

Whole-heart RNA-seq

RNA-seq of hearts from P16 littermate control and knockout was as described previously (Wang et al. 2017). The data were deposited in the GEO database (GSE88761).

RNA and protein analysis

qRT–PCR, luciferase reporter assay, and immunohistochemistry were performed as described previously (Pei et al. 2006, 2011, 2015; Wang et al. 2017; Zhao et al. 2018). We isolated total RNA from mouse tissues or HL-1 cells using RNAzol RT (Molecular Research Center) following the manufacturer's instructions. The antibodies used were PECAM1 (Santa Cruz Biotechnology, sc-46694; and BD Biosciences, 550274), NPR3 (Santa Cruz Biotechnology, sc-515449), TNNI3 (Abcam, ab47003), Ki67 (Vector Laboratories, vp-RM04), and GDF15 (Abcam, ab189358). Sirius Red staining was performed using the Picrosirius Red stain kit (Polysciences). Adenovirus construction and pericardial injection were performed as described previously (Pei et al. 2006, 2011; Wang et al. 2017).

Statistical analysis

Fishers’ exact test, χ2 test, correlation test, and Student's t-test were performed in R to determine the statistical significance, with P < 0.05 deemed as statistically significant.

See the Supplemental Material for additional Materials and Methods.