Transcriptional heterogeneity of fibroblasts is a hallmark of the aging heart

In general, aging increases the incidence of cardiovascular disease and significantly alters cardiac structure and function (1). Intervention with cardioprotective substances has been proposed to extend life span (2). Cardiovascular aging is associated with increased vascular and ventricular stiffness, interstitial cardiac fibrosis, augmented left ventricular wall thickness, and reduced neovascularization abilities after ischemia (1, 3). Fibrotic remodeling of the aging ventricle plays an important role in the pathogenesis of diastolic heart failure, and fibrosis of the conduction system contributes to the development of arrhythmias and conduction abnormalities (4, 5). On a cellular level, replicative senescence caused by shortening of telomeres, increased oxidative stress and DNA damage, deregulation of genes and proteins, impaired cell-cell communication, and an altered systemic and local environment are hallmarks of aging (6). In the heart, aging is associated with increased fibroblast activation, impaired angiogenesis, and increased cardiomyocyte hypertrophy, together contributing to the above-described cardiac phenotypes (7). In addition, relative numbers and proportions of cardiac leukocytes are altered: old hearts have proportionally more monocyte-derived cardiac macrophages and an increased population of granulocytes (8). However, we still know surprisingly little about the fundamental effects of aging on cellular heterogeneity of cardiac parenchymal, mural, and vascular cells.

Here we describe what may be one of the first mammalian transcriptional cell atlases of the aging heart. These data allow understanding of both the healthy state and the causal effects of intrinsic cardiac aging.

Single-nucleus RNA-sequencing reveals cellular heterogeneity of cardiac cells in aging. To comprehensively decipher the expected cellular responses to intrinsic cardiac aging, we applied microdroplet-based single-nucleus RNA-sequencing (9) of cross-sections of snap-frozen heart samples from 3 syngeneic young male mice (12 weeks old) and 3 aged male mice (18 months old). In total, we analyzed 14,827 nuclei from young hearts and 12,981 nuclei from old hearts (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.131092DS1). Using t-distributed stochastic neighbor embedding (tSNE) (10), global dimension reduction was constructed from all samples to visualize clusters that were defined by cell-specific gene markers (Figure 1A and Supplemental Table 2). Alignment of samples indicated high reproducibility across samples (Supplemental Figure 1). Most of the cells were in G₁ phase, and no influence of aging on cell cycle activity was observed (Supplemental Figure 2).

Figure 1

Cell type–specific composition of the heart in young and old mice. (A) tSNE representation of all cells and cell clusters from all 6 samples. Clusters were annotated according to their detected gene markers. Details about each cell population can be found in Supplemental Figure 3 and Supplemental Table 2. Most cell clusters could be separated with high confidence; however, in this initial analysis, fibroblasts remained difficult to distinguish from one another. (B) A total of 128 unique differentially expressed genes (DEGs) (adjusted P < 0.1) were found between young and old samples among all detected clusters. Outer circle represents upregulated genes in old samples, and inner circle represents the downregulated genes in old. (C) GO enrichment comparison (hypergeometric test) of the DEGs between young and old samples in the cell populations with at least 1 significant result (adjusted P < 0.1). Up- and downregulated genes were analyzed together. Subpopulations were analyzed together. (D) The DEGs were grouped into coexpressed networks and represented as different colors; these networks were functionally annotated according to their genes. These genes were spatially organized in a Venn diagram for easy access of same DEGs in multiple cell types.

Unsupervised clustering revealed 15 distinct gene expression patterns (Figure 1A and Supplemental Figure 3). Using cell type–specific gene markers (Supplemental Table 2) and published mouse single-cell gene expression data (11, 12), 7 major cell types could be annotated, including fibroblasts (A, B), cardiomyocytes (A, B, C), endothelial cells (A, B, C), immune cells (A, B, C), pericytes, epicardial cells, and adipocytes (Figure 1A and Supplemental Figure 3). In particular, for fibroblasts, the unsupervised clustering revealed 2 main clusters, fibroblast A (79.42%) and fibroblast B (20.58%). Separation of these 2 clusters was not significant (Supplemental Figure 3B), and gene markers were very similar (Supplemental Table 2); moreover, these 2 clusters were almost equally populated by young and old cells. Analysis of the cell numbers in clusters of other cell types than fibroblasts showed in part trends for changes during aging (Supplemental Figure 4) but did not reveal statistically significant differences.

In general, 128 differentially expressed nonredundant genes (DEGs) were found between young and aged hearts (Figure 1B and Supplemental Table 3). Considering the DEGs in all cell clusters, 107 genes showed significantly increased expression (adjusted P < 0.1), and 21 genes showed significantly decreased expression (adjusted P < 0.1) in aged versus young hearts (Supplemental Table 3). Interestingly, aging predominantly affected gene expression patterns in fibroblasts (Figure 1B). Several highly differentially expressed genes could be confirmed by quantitative reverse transcription PCR of isolated cardiac fibroblasts (Supplemental Figure 5).

Gene Ontology (GO) analysis of DEGs revealed a cell type–specific enrichment of genes associated with various pathways, such as angiogenesis, chemotaxis/migration, inflammation/immune response, and cell/matrix association (Figure 1C). Only a few coexpression networks and significantly regulated genes were shared between the main cell types. Among them, the expression of the components of the complement system were commonly augmented in all cell types (Figure 1D, Supplemental Table 4), which is consistent with the finding of a general cardiac aging-promoting effect of the complement system (13).

Single-nucleus RNA-sequencing identifies specific fibroblast subpopulations involved in cardiac aging. Because our data suggest that aging has the most profound impact on cardiac fibroblasts (Figure 1, B and D), we focused our attention on these cells. To gain insights into age-associated fibroblast populations, we applied subclustering techniques to sort and group cells using the 85 unique genes that were differentially expressed in fibroblasts during aging. Subclustering identified 12 fibroblast subpopulations (Figure 2, A and B), of which subclusters 1, 5, 10, and 11 were mostly populated by young cells, whereas subclusters 2, 3, 4, 7, 8, and 12 were mostly populated by old cells (Figure 2C). All gene markers for the subclusters are listed in Supplemental Table 5.

Figure 2

Aging fibroblast subclusters. (A) Old and young fibroblasts were clustered using DEGs as anchors. (B) Subclusters of fibroblasts were found using Louvain algorithm implemented in Monocle2 (see Methods). (C) Differences in cell numbers between old and young in all subclusters were calculated by comparing the normalized means and the significance by Fisher’s exact test (*P < 0.05; **P < 0.01). (D) GO enrichment (hypergeometric test with adjusted P value cutoff of 0.1) comparison between all fibroblast subclusters (if any).

The different age-dependent subclusters showed profound alterations in gene pathways (Figure 2D and Supplemental Table 6). GO terms that were enriched in subclusters of young fibroblasts were connected to regulation of muscle system process and generation of metabolites (subcluster 11) (Figure 2D). Angiogenesis was in general found to be enriched in all subclusters that were mainly populated by old fibroblasts. Moreover, subclusters 2, 4, and 8, also mainly populated by old cells, associated with genes involved in regulation of endothelial cell proliferation and migration. Genes associated with osteoblasts were also prominently augmented in subclusters 2 and 4. Immune response genes were mainly found in subcluster 3 (Figure 2D and Supplemental Table 6), which was predominant for cells from old mice. The increase in proinflammatory phenotypes resembles the response of fibroblasts during cardiac injury or remodeling (4, 5). The age-associated phenotype of subcluster 2 showing an induction of epithelial cell proliferation is consistent with an expansion of fibroblasts and increased interstitial fibrosis in the aging heart (ref. 4; Supplemental Figure 6). However, it was unknown that cardiac fibroblasts elicit cell-specific individual changes in vascular/angiogenesis profiles during aging (as observed in subclusters 2, 3, 4, 7, and 8) or may acquire more cells with osteogenic phenotypes (as in subclusters 2 and 4).

Single-nucleus RNA-sequencing reveals dynamic states of fibroblasts and entropy increase during aging. Next, we applied the DDRTree algorithm from Monocle2 to sort and plot cells in a linear order. The linear representation of the cells shows the dynamic states of fibroblasts during the aging process (Figure 3, A and B, and Supplemental Figure 7). Here, states of the right end (states a, b, and d) were dominantly populated by young fibroblasts, which showed increased expression of genes required for heart fibroblast functions, such as extracellular structure organization. States of the left end (states h, j, and l) were predominantly populated by aged fibroblasts, which were enriched with genes required for cell proliferation (j), inflammatory response (h, j, l), angiogenesis (j, l), and osteoblast differentiation (j, l) (Figure 3, A and C, Supplemental Table 7, and Supplemental Table 8).

Figure 3

Transcriptional activity in aging cardiac fibroblasts. (A) Fibroblasts were sorted using DDRTree algorithm and using the DEGs between old and young to define different cell development states during aging. (B) Projection of subclusters onto the different cell states. Same colors from Figure 2B were used. (C) Differences in cell numbers between old and young of all fibroblast states were calculated by comparing the normalized means. Significance was calculated by Fisher’s exact test (**P < 0.01). (D) Transcriptional activity was estimated by measuring the ratio between unspliced and spliced mRNAs. This so-called RNA velocity is represented by high-dimensional vectors; the longer the arrow in the plot, the higher the transcriptional activity as seen in the extremities of states h, j, and l plot containing mostly old cells.

Interestingly, analyzing the entire single-nucleus mRNA expression data of cardiac fibroblasts, we observed high numbers of unspliced RNAs (increased transcriptional activity or “RNA velocity”; ref. 14) in the fibroblast states j, l, and h (Figure 3D); these most dynamic states were particularly populated by old cells (Figure 3, A–D).

Furthermore, we observed a consistent age-related increase of entropy in all fibroblasts when comparing young and old fibroblasts (Figure 4A). Subpopulations that were mainly populated by old cells showed higher entropy than subpopulations that were rather populated by young cells (Figure 4B). Of note, genes related to the cell-specific functions of fibroblasts, such as extracellular matrix organization, showed increased entropy in aged versus young fibroblasts (Figure 4C). These results indicate a higher transcriptional heterogeneity of aged heart fibroblasts, which could lead to a decline of efficient, synchronous expression of fibroblast-specific genes, eventually resulting in decreased molecular network stability.

Figure 4

Single-cell entropy. (A) Single-cell entropy observed in old and young fibroblasts. (B) Single-cell entropy (transcriptional “disorder”) detected in fibroblast subclusters. Highest entropy was found in subclusters 2, 4, 7, and 8, which were mostly populated by aged fibroblast cells. (C) Single-cell entropy observed in overall old and young fibroblasts for selected gene sets. Mann-Whitney-Wilcoxon test was applied in the pairwise comparisons (young versus old) to calculate statistical significance between young and old cell entropy means (*P < 2.2e-16). The box plots depict the minimum and maximum values (whiskers), the upper and lower quartiles, and the median. The length of box represents the interquartile range.

Aging alters fibroblast/endothelium interactions. To analyze potential functional consequences of the observed transcriptional changes, we first assessed the differentially expressed angiogenesis-related genes found in fibroblast cell clusters 2, 3, 4, 7, and 8 (Figure 2, B and D). Interestingly, many of these genes are known from endothelial-mesenchymal transition (EndMT), which refers to the process leading to the phenotypic change of endothelial cells into mesenchymal cells (ref. 15; Figure 5A). Therefore, we speculated that endothelial cell–derived fibroblasts might have accumulated during aging. However, no significant increase of endothelial markers was found in old fibroblasts, suggesting that the EndMT-associated regulated genes (e.g., Klf4, Tgfbr2, Sulf1, Loxl2; Figure 5B) represent the activation of a transcriptional transition program in cardiac fibroblasts but are not related to the transition of the endothelial lineage.

Figure 5

Fibroblast and endothelium interaction. (A) Genes associated with endothelial-mesenchymal transition (EndMT) and epithelial-mesenchymal transition (EMT) expressed in fibroblast states l, j, and h. (B) Expression of Klf4, Loxl2, Sulf1, and Tgfbr2 in fibroblast subclusters. (C) Interaction network of each annotated cluster showing the total outgoing weights (secreted/paracrine/autocrine signals) and total incoming weights (ligand binding/receptor activation) signals among the clusters. (D) Ligand/receptor interaction analysis. Connections found between fibroblasts and endothelial cells are shown for young and old mouse hearts.

Some of the observed angiogenesis-related genes associated with aging are also known to encode for extracellular proteins, suggesting that a paracrine crosstalk of cardiac fibroblasts with cardiac endothelial cells may be affected during aging. Indeed, weighted analysis of the ligand/receptor interactions indicated that fibroblasts were the most “outbound” cells, whereas endothelial cells were the most “receiving” cells in the aging heart (Figure 5, C and D). To address whether the paracrine activity of fibroblasts is altered in aging, we generated conditioned medium from aged and young mouse cardiac fibroblasts and determined the angiogenesis regulatory activity in vitro (Figure 6A). Conditioned medium derived from aged fibroblasts showed a reduced angiogenic property compared with medium of young heart-derived fibroblasts, demonstrating an impaired crosstalk during aging (Figure 6B). Consistently, capillary density was profoundly reduced in fibrotic regions of aged hearts (Figure 6C and Supplemental Figure 8A).

Figure 6

Age-dependent changes of fibroblast function on endothelial cells. (A and B) Experimental outline and tube formation assay of human umbilical vein endothelial cells (HUVECs) that were cultured in conditioned medium received from young (12 weeks) and aged (20 months) cardiac mouse fibroblasts. Accumulated tube length was measured in 5 randomly chosen microscopic fields with a computer-assisted microscope using Axiovision 4.5 (Zeiss) (scale bar: 100 μm; n = 4). (C) Capillary density of random healthy areas of young (12 weeks; n = 6) and old (n = 8) heart sections versus capillary density of fibrotic (n = 8) areas of aged hearts (20 months). Shown is the quantification of the Isolectin B4–positive area versus the total area (%). (D) Expression of Serpine1, Serpine2, and Efemp1 in fibroblasts displayed in log scale in trajectory plots. (E) Secretion of SerpinE1 in isolated young and aged cardiac fibroblasts cultured for 24 hours in serum-free medium. Data were derived using a mouse angiogenesis array (R&D Systems) of culture supernatants of isolated cardiac fibroblasts (n = 4). (F) Tube formation assay of HUVECs that were cultured in conditioned medium received from young (12 weeks) and aged (20 months) cardiac mouse fibroblasts. The aged phenotype was rescued by supplementing the aged fibroblast medium with 4 μg of an anti-Serpin antibody. Accumulated tube length was measured in 5 randomly chosen microscopic fields with a computer-assisted microscope using Axiovision 4.5 (Zeiss) (scale bar: 200 μm; n = 4). (G) RNA expression of IL-6 in HUVECs cultured in conditioned medium received from young (12 weeks) and aged (18 months) cardiac mouse fibroblasts and in the presence and absence of 10 ng/mL TNF-α for 24 hours (n = 5). (H) Quantification of autophagosomes in HUVECs cultured with conditioned medium received from young (12 weeks) and aged (18 months) cardiac mouse fibroblasts. Autophagosomes were detected by immunostainings against p62 and LC3, whereby autophagosomes were considered p62/LC3–double positive (n = 5). Data are shown as mean ± SEM (B, C, E, F, G, and H). After passing Gaussian distribution, statistical analysis was performed using the unpaired, 2-sided t test (C, E, G, and H). For comparisons of more than 2 groups, multiple-group 1-way ANOVA with a post hoc Bonferroni’s test was used (B and F). *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Aged fibroblasts expressed various candidates that might mediate antiangiogenic effects, including Efemp1 (also known as Fibulin 3), a gene that is known to inhibit tumor angiogenesis (16), and several members of the antiangiogenic serpins such as Serpine1 (also known as PAI1) and Serpine2 (alias Nexin) (ref. 17; Figure 6D). Using a screen for angiogenesis-related proteins, we detected increased levels of Serpine1 in the medium of aged fibroblasts (Figure 6E). Next, we tested the effect of recombinant Serpine1 and Serpine2, which both reduced tube formation (Supplemental Figure 8, B and C). To determine their contribution to the antiangiogenic activity of aged fibroblast supernatants on endothelial cells, we used neutralizing antibodies that block both serpins (Supplemental Figure 8, B and C). Inhibition of serpins significantly rescued the impaired angiogenesis induced by conditioned medium from aged fibroblasts (Figure 6F), whereas Efemp1 antibodies had no effect on endothelial cells (Supplemental Figure 8D). These data suggest that increased expression of serpins in fibroblasts during aging mediates an antiangiogenic effect on endothelial cells.

Effect of fibroblasts on age-associated signatures of endothelial cells. To determine whether aged heart-derived fibroblasts may further promote an aging phenotype of endothelial cells, we tested the effect of conditioned medium on hallmarks of aged endothelial cells. Conditioned medium of aged heart-derived fibroblasts did not affect proliferation of endothelial cells (Supplemental Figure 8E), suggesting that it does not influence replicative senescence. However, aged fibroblast-conditioned medium enhanced the proinflammatory activation of endothelial cells (Figure 6G) and reduced endothelial cell autophagy (Figure 6H), which both are typical features of impaired endothelial cells (18, 19). These data suggest that aged heart fibroblasts impair endothelial cell proangiogenic functions and promote age-associated changes in endothelial cells.

Aging fibroblasts adopt osteogenic fates. Finally, we addressed the relevance of the observed augmented and more entropic expression of osteogenic genes in fibroblasts of subclusters 2 and 4 (Figure 2D, Supplemental Table 6, and Figure 7A). Cardiac fibroblasts were shown to adopt osteogenic fates in conditions of pathological heart remodeling (20); however, whether such transitions occur during aging remained unclear. Therefore, we further explored these findings by analyzing the expression of genes associated with osteogenesis. Various genes involved in osteoblast differentiation and functions (e.g., Ddr2, Runx2, Gpm6b, Junb, Cebpb) were all expressed at higher levels in subcluster 2, which was populated by old fibroblasts (Figure 7A and Supplemental Table 5). Some of these genes, such as Runx2, were also significantly elevated in the total fibroblast population (Figure 7B). Histological assessment confirmed the occurrence of a subset of osteocalcin-expressing fibroblasts in the aging heart (Figure 7C), which were mainly located in the epicardial layer (Figure 7, D and E). Consistently, calcification was increased in the aging heart, particularly in the epicardium (Supplemental Figure 9A). To determine whether epicardial cells may acquire an osteogenic state, we determined the expression of the osteoblast genes in the epicardial cell cluster (Supplemental Figure 9, B and C). The epicardial cell cluster is characterized by high expression of epicardial markers Wt1 and Upk3b but only weakly expresses fibroblast marker genes, such as Pdgfra (Supplemental Figure 9B). However, this cluster did not show any changes in the expression of osteogenic genes, such as Cebpb or Runx2 (Supplemental Figure 9C), suggesting that the acquirement of the osteogenic state may preferentially occur in Pdgfra⁺ fibroblastic populations within the epicardial layer.

Figure 7

Aging alters fibroblasts to adopt osteogenic fates. (A) Expression of Cebpb and Runx2 in fibroblasts displayed in log scale in subcluster analysis plot. Young and old enriched subclusters are separated by lines and as displayed in Figure 2, A and B. (B) Expression of Runx2 in young and aged cardiac fibroblasts. Data were derived from bulk RNA-sequencing of isolated cardiac fibroblasts (4 young vs. 3 old samples). FPKM, fragments per kilobase per million mapped reads. (C–E) Histological analysis of heart sections derived from young (12 weeks old) and old (>18 months old) mice. Immunostaining was performed against osteocalcin (OCN, green) as an osteogenic marker and PDGFR-α (gray) as a fibroblast marker. Nuclei were stained with Hoechst 33342 (blue) (n = 8; 4–6 images/heart). Quantification of OCN⁺PDGFR-α⁺ double-positive cells vs. PDGFR-α⁺ single-positive cells is shown in C. Representative examples are shown in D (scale bar: 10 μm) and E (scale bar: 100 μm). Images were monitored using the Leica TCS SP8 confocal microscope. Data are shown as mean ± SEM (B and C). After passing Gaussian distribution, statistical analysis was performed using the unpaired, 2-sided t test. *P < 0.05; ****P < 0.0001.

Here we analyzed transcriptional changes of major cell types of the aging mouse heart at single-cell resolution. Of note, fibroblasts showed most profound heterogeneity and the highest number of DEGs. In particular, in cardiac fibroblasts, we observed more fuzzy and diverse transcriptional profiles in old versus young heart cells, as evidenced by increased entropy of cell type–specific, functionally relevant genes and pathways, indicating quantitative loss of systemic function of fibroblasts during cardiac aging. These results suggest a general age-associated change in cellular heterogeneity of cardiac fibroblasts, which is in line with previous studies showing an age-dependent increase in cell-cell transcriptional variability after immune stimulation of CD4⁺ T cells in mice (21). Possibly, the increased RNA velocity of fibroblasts may relate to resolution of cell-specific boundaries associated with dedifferentiation or acquirement of interim cellular phenotypes.

Fibroblasts are well known to contribute to interstitial fibrosis, which is a pathological process that increases the stiffness of the heart and contributes to age-associated decline in diastolic heart function. Interstitial fibrosis is associated with deregulation of resident mesenchymal fibroblasts, which acquire a proinflammatory state (22). This finding is recapitulated in our trajectory analysis, demonstrating immune/inflammatory activation in the fibroblasts, in particular in states j and l, which are preferentially populated by cells from old heart-derived fibroblasts.

Various previous studies suggest that fate changes between endothelial cells and fibroblasts can occur in both directions under pathological conditions. Both EndMTs and mesenchymal-endothelial transitions were described after pathological activation of the heart, e.g., after infarction. The incidence of EndMT in aging seems largely unknown, but senescent endothelial cells were shown to acquire mesenchymal properties in vitro (23). However, we could not find evidence for an increase of cells coexpressing mesenchymal and endothelial markers, which is compatible with recent lineage tracing studies, suggesting also that the extent to which endothelial cells contribute to heart fibroblasts after myocardial infarction is rather limited (24).

Our comprehensive single-cell transcriptome data and weighted analysis of the ligand/receptor interactions suggest interaction of fibroblasts and “receiving” endothelial cells in the aging heart. Single-nucleus RNA-sequencing data indicated deregulation of angiogenesis-related genes with known paracrine functions in the aging fibroblast cell populations. In line with this, our functional tests revealed a disturbed crosstalk of aged cardiac fibroblasts with endothelial cells during aging. Fibroblasts may contribute to the change of the microenvironment in the aged heart to control the functional response of endothelial cells. Specifically, we identified fibroblast-derived serpins as mediators of endothelial dysfunction. Serpine1 expression is well known to be increased in senescent cells in the aging heart (25). A null mutation of SERPINE1 in humans protects against aging (26). In combination with its antiangiogenic activity (17), the induction of this pleiotropic gene in aging fibroblasts may contribute to the reduced capillary density observed in the fibrotic regions of the aging hearts. The impaired microcirculation observed in the aging heart may further aggravate fibrosis, cardiac hypertrophy, and cardiac dysfunction. In addition, we showed that conditioned medium increased the proinflammatory activation of endothelial cells and reduced endothelial cell autophagy. While the mechanism by which aged fibroblasts may interfere with endothelial cell autophagy is unclear, one may speculate that Serpine1, which was shown to exhibit a proinflammatory activity (27), may also affect endothelial inflammatory responses. Enhanced vascular inflammation may lead to tissue invasion of inflammatory bone marrow–derived cells, which may replace tissue-resident macrophages and further enhance the aging process. Interestingly, conditioned medium from aged fibroblasts had no significant effect on cellular proliferation, suggesting that aged fibroblasts do not contribute to endothelial cell replicative senescence. Together, our data suggest that fibroblasts of the aging heart induce some but not all the signatures of aging endothelial cells. In addition to serpins, many other cytokines and extracellular matrix genes are dysregulated in the aging fibroblasts, many of which are predicted to interact with endothelial cells. Therefore, it is likely that the interaction of aging fibroblasts with endothelial cells in vivo is more complex than could be mimicked by mechanistic experimental intervention.

In addition, we found that old fibroblasts of the heart acquired osteogenic traits. Interestingly, a very recent study showed that aged skin fibroblasts appear to gain adipogenic traits (28), suggesting a tissue-specific impact of aging on fibroblast identity. Aging is well known to augment vascular calcification, but a fate change of fibroblast to osteogenic cells in the heart has only recently been described under pathological stress conditions (20). Our data suggest that natural aging is also associated with such a change in fibroblast fates, leading to an increasing, but low, number of osteocalcin gene–expressing cells in the heart. Previous published single-cell data sets of young mouse hearts (11, 29) showed mostly cells with gene expression profiles resembling here-shown young fibroblast subclusters (Supplemental Figure 10). Very few cells from these public data sets have transcriptional profiles related to fibroblast subclusters 2, 4, or 8 (enriched with genes involved in epithelial cell proliferation and osteoblasts). Histological staining suggested that osteocalcin-expressing cells were located in the epicardial layer, in which calcification in particular was observed. Our in silico analysis of the epicardial cell cluster suggested that aged heart epicardial cells do not show an increase in osteoblast genes, but we cannot rule out that a subset of Pdgfra-expressing epicardial cells of the aging heart may cluster within the fibroblast population. The consequences of these fate changes on cardiac function need to be further explored.

In this study we focused on the discovery of cell type–specific effects of aging under well-defined conditions. Our single-nucleus transcriptome data revealed systemic responses of cardiac cells to aging, in particular with respect to the heterogeneity and specific age-related changes of fibroblast subpopulations (Supplemental Figure 11). Of note, beyond shedding light on mechanistic processes of aging with unprecedented resolution, transcriptome analyses of individual cells can be leveraged to discover entry points of physiological disorders resulting from the effects of intrinsic aging.

Laboratory animals. Isogenic male C57BL/6JRj inbred mice were obtained from Janvier (Le Genest Saint-Isle, France). Homozygosity of these inbred mice was controlled by Janvier using exome sequencing (https://www.janvier-labs.com/461.html). To obtain hearts, mice were sacrificed by cervical dislocation under anesthesia with isoflurane. The obtained hearts were either snap-frozen with liquid nitrogen for nuclei isolation or freshly used for cell isolation.

Nuclei isolation. Whole mouse hearts (derived from three 12-week-old and from three 18-month-old mice) were dissociated in liquid nitrogen using a mortar. The sample was then transferred into a Dounce homogenizer. For isolation of nuclei we strictly followed the experimental steps that were detailed previously (30). Samples were homogenized and lysed to release the nuclei, followed by FACS (Aria III, BD Genomics) of DAPI-stained high-quality nuclei to sort out any residual cell debris. Integrity and purity of the nuclei were inspected visually using a Leica DMi8 inverted microscope. Yield of nuclei was quantified by using a Neubauer chamber.

Single-nucleus cDNA library preparation and Illumina short-read sequencing. Generation of libraries for single-cell mRNA-sequencing was done using the Chromium Single Cell 3′ protocol (10x Genomics). Briefly, we applied droplet-based nuclei encapsulation, lysis, as well as individual transcript (UMI) and nuclei barcoding for generation of 3′ end-counting mRNA-sequencing libraries, using the Single Cell 3′ v2 Reagent Kits from 10x Genomics. The libraries were sequenced using a HiSeq 4000 sequencer (Illumina). Bulk transcriptome data of fibroblasts generated by applying conventional TruSeq library preparation (Illumina) and Illumina sequencing were used to validate single-nucleus sequencing results.

Quality control and read mapping. The Cell Ranger suite version 2.1.0 was used to perform barcode processing and single-cell gene UMI counting (http://software.10xgenomics.com/single-cell/overview/welcome).

Reads were aligned to GRCm38 reference genome, and gene annotation was extracted from gencode vM12. Intronic mapped reads were assumed to be part of pre-mRNAs; therefore, they were also assigned to the respective gene. Because Cell Ranger counts only those reads that map to exons, we created a modified gene transfer format file, which denotes each gene transcript locus as 1 exon and rebuilt a Cell Ranger reference (cellranger mkref). This step was necessary because of a high number of pre-mRNA in nuclei. Final results of the Cell Ranger analysis contain the count values of UMIs assigned to each gene in each of the cells for each respective sample using all mapped reads.

Cell cycle analysis. Cell cycle states were established using cyclone implemented in scran package (v1.6.6) and the mouse model available in the package (mouse_cycle_markers.rds). Cells were classified into G₁, S, and G₂M phase, based on their normalized gene count (31).

Data normalization, clustering, and cell population definition. Cell Ranger raw result tables were merged and analyzed using Seurat suite version 2.2.0 (32). Only genes expressed in at least 40 different nuclei and nuclei with a minimum of 200 genes were generally kept for analysis. For some specific analyses below, we partly applied more restricted values. Gene expression (in UMI) was log normalized and scaled. Variable genes for each of the 6 samples were detected by the Seurat “FindVariableGenes” function using default options. Clustering was done using only variable genes. To obtain a global set of variable genes, we combined the set of variable genes from all 6 samples using the union of the top 1000 variable genes in each one. To align all 6 samples and remove potential batch effects derived from sample preparation and sequencing, we first ran a canonical correlation analysis (CCA) using a diagonal implementation of CCA on each sample. Then we aligned the subspaces across all selected variables on this CCA.

To explore transcriptional heterogeneity and to undertake initial cell clustering, we reduced dimensionality using the aligned CCAs. We used CCA loadings as input for a graph-based approach to cluster cells by cell type and as input for tSNE for reduction to 2 dimensions for visualization purposes. For the overall data set, we selected the first 20 dimensions that explained more variability than expected by chance using a permutation-based test as implemented in Seurat. For further clustering of cells within primary clusters (subclustering), we selected variable numbers of CCAs for dimensionality reduction using either a permutation-based test or heuristic methods implemented in Seurat. Clusters were identified using the function “FindClusters” from Seurat using default parameters. The robustness of the clusters was calculated using the function “AssessNodes” from Seurat. For each cluster, a phylogenetic tree based on the distance matrix in gene expression space was computed. Next, Seurat computed an out-of-bag error for a random forest classifier trained on each internal node split of the tree).

Cluster annotation to heart cell types. To identify specific markers for each cluster and consequently annotate them with the respective heart cell type, we applied Seurat “FindAllMarkers” function using the AUC-based scoring classifier with default parameters. We used the default Wilcoxon’s rank-sum test for the specific cell type analysis.

The assignment of each cluster to a specific cell type was done manually by making use of published gene expression data and by single-cell profile projection using scmap. We projected single heart cells from the Tabula Muris consortium (12) set onto our single nuclei to identify the similar cell types and show the reproducibility from both approaches. All default parameters were used; the numbers of selected features were defined as the top 40.

Differential expression using DESeq2. To identify DEGs between the old and the young hearts, we used DESeq2 (v. 1.20) to explore the statistical power of the replicates (33). Differential expression analysis was conducted separately for each cluster. Data were filtered for quality; for each sub cell type, we allowed only cells/nuclei with more than 1000 UMIs (for cardiomyocytes we used 2000 as minimum because of the high number of nuclei from this type), 500 detected genes, and maximum 0.05% of mitochondrial reads. Only clusters with more than 500 cells after filtering were analyzed. Clusters were summarized by total counts for each gene to generate an in silico bulk RNA-sequencing data set. This way we could assess fold changes and adjusted P values for DEGs between young and old samples, containing 3 replicates for each subtype. The test for significance was computed by the nbinomLRT DESeq2 function, and the significantly differentially expressed genes were selected using independent hypothesis weighting (34).

Gene enrichment analysis. We used ClusterProfile (R package; v. 3.8) and ConsensusPathDB (http://cpdb.molgen.mpg.de/CPDB) for pathway analysis of age-associated transcripts. In ClusterProfile, all gene names were first converted into Entrez IDs and then enriched for both Kyoto Encyclopedia of Genes and Genomes and GO terms. Benjamini-Hochberg FDR was used for multiple-testing corrections. The significance threshold of P < 0.05 after correction for multiple testing was applied. We further used the online tool Metascape (http://metascape.org/gp/index.html#/main/step1) for GO analysis of genes enriched and differentially expressed in fibroblast populations A and B.

Coexpression networks. WGNCA (35) was used to build coexpression networks from DEG profiles. Each cluster was analyzed independently, as were up- (higher expression on old hearts) and downregulated genes. Cells were extra filtered for quality as described in Differential expression using DESeq2. The soft threshold power parameter was selected when “scale-free topology r-squared” reached a value higher than 0.7. Minimum module size was set to 5.

Subclustering based on DEGs. DEGs between old and young fibroblasts were used for subclustering and detecting cell populations of fibroblasts based on the signal produced by these genes. We used Monocle2 sorting cell functions and set up the DEGs as ordering genes for reducing dimensions. Seurat objects were converted into Monocle (v2.4) (36) objects using the “importCDS” function. Two dimension reduction algorithms were used for graphical representation of the data, the tSNE and DDRTree. For tSNE we used the following parameters: num_dim = 15, max_components = 2, and norm_method = log. The clustering was performed using the following parameters: rho_threshold = 2, and delta_threshold = 5. For DDRTree we used this as the parameter: auto_param_selection = T.

The DDRTree graph representation helped us visualize the progress of possible cell transitions within the cell types using aging as the sorting factor. We plotted the subclusters split by age and states split by subcluster and age.

RNA velocity. We ran velocyto.py (14) annotator for each mapped BAM file using the default parameters for 10x Genomics technology and gencode vM12 gtf file for intron-exon annotation. The resulting loom object for each sample was loaded and processed in R using the velocyto.R (v. 0.17) package. We used the embedding from Monocle DDRTree representation for cell-cell distance calculation and final velocity plots. The estimation of RNA velocity was done with the grouping of 10 cells; the other parameters were set to default.

Entropy calculation. We used SLICE (v0.99.0) (37) package to calculate entropy values of groups of cells/nuclei based on a set of genes. The Seurat objects were sliced according to this set of genes and to the group of cells/nuclei to be analyzed. Entropy was calculated by using the getEntropy function from SLICE package applying a bootstrap calculation of size 1000, 100 iterations, and random seed “201602.”

Ligand/receptor interactions. Ligand/receptor interaction analysis was performed according to the approach described previously (29). Briefly, a weighted directed graph with 4 layers of nodes was built linking source cell types (layer 1), defined by expression of a ligand (layer 2), to target cell types (layer 4) expressing a corresponding receptor (layer 3), after reference to a curated map of ligand/receptor pairs (38). Source/ligand and receptor/target edges were weighted according to expression fold-change in ligands and receptors, respectively. Ligand/receptor edges were further weighted by mouse-specific association scores from the STRING database (39). Permutation testing (100,000 permutations) of randomized network connections was applied to determine significant source-target network connections following Benjamini-Hochberg multiple-testing correction (adjusted P < 0.01).

Cardiac fibroblast cell isolation. Cardiac fibroblasts were isolated from male old (20 months) versus young (12 weeks) mice. After sacrificing the mice, hearts were harvested, cut into small pieces, and washed with Hanks’ buffered saline solution (+Ca²⁺/+Mg²⁺). Tissue dissociation was performed in 2.5 mL of a commercial enzyme mix (Neonatal Heart Dissociation Kit, mouse and rat, Miltenyi Biotec GmbH). To dissociate the solid heart tissue, gentleMACS Dissociator (Miltenyi Biotec GmbH) with the preprogrammed program m_neoheart_01_01 was used after 25 minutes, 3 times 15 minutes of digestion at 37°C. The cellular components of the digested heart suspension were pelleted by centrifugation (80 g, 5 minutes, 4°C) and then resuspended in 3 mL DMEM High Glucose GlutaMAX (Invitrogen, Thermo Fisher Scientific) containing 10% FCS (Invitrogen, Thermo Fisher Scientific) and penicillin/streptomycin (Roche). Fibroblasts were isolated by plating the cell suspension for 2 hours at 37°C and 5% CO₂ in a humidified atmosphere. After plating, the nonfibroblast cells were removed by 3 washing steps using PBS (Invitrogen, Thermo Fisher Scientific).

RNA isolation and quantification. Total RNA was purified from cells using miRNeasy kits (Qiagen), combined with on-column DNase digestion (DNase Set, Qiagen) according to the manufacturer’s instructions. The RNA concentration was determined by measuring absorption at 260 nm and 280 nm with the NanoDrop ND 2000-spectrophotometer (PeqLab).

Bulk RNA-sequencing. RNA was isolated from isolated cardiac fibroblasts using the miRNeasy Micro Kit (Qiagen) combined with on-column DNAse digestion (DNase-Free DNase Set, Qiagen) to avoid contamination by genomic DNA. RNA and library preparation integrity were verified with a BioAnalyzer 2100 (Agilent) or LabChip Gx Touch 24 (PerkinElmer). Two hundred nanograms of total RNA was used as input for TruSeq Stranded Total RNA Library preparation following the low sample protocol (Illumina). Sequencing was performed on the NextSeq500 instrument (Illumina) using v2 chemistry, resulting in a minimum of 30 million reads per library with 1 × 75 bp single-end setup. Demultiplexing and adapter removal were performed by bcl2fastq software (Illumina), and resulting raw reads were assessed for quality, adapter content, and duplication rates with FastQC. Fastx_trimmer was used to trim the first 5 bp. Trimmed and filtered reads were aligned versus the Ensembl mouse genome version mm10 (GRCm38) using STAR 2.4.0a with the standard parameters (40). The number of reads aligning to genes was counted and compared with Ciffdiff version 2.2.1 (41). The Ensembl annotation was enriched with UniProt data (release 06.06.2014) based on Ensembl gene identifiers (Activities at the Universal Protein Resource [UniProt]).

cDNA synthesis. To quantify mRNA expression by qPCR, 100 ng to 1 μg of total RNA was reverse-transcribed using the reverse transcriptase MulV (Life Technologies, Thermo Fisher Scientific) as previously described (42). The cDNA synthesis was performed in a reaction volume of 20 μL using 1× PCR Buffer II (10×) with magnesium (Thermo Fisher Scientific), 5 mM of MgCl₂ (Applied Biosystems), 0.5 μg of random hexamer primer (Thermo Fisher Scientific), 0.5 mM of dNTP mix (Fermentas), 20 units of RNAse inhibitor (Thermo Fisher Scientific), and 50 units of MulV reverse transcriptase. The reaction was incubated at 43°C for 75 minutes followed by 5 minutes of heat inactivation at 95°C. The synthesized cDNA was diluted with RNAse-/DNAse-free water (Invitrogen, Thermo Fisher Scientific) to prevent any interference during quantitative PCR caused by the ion strength present in the cDNA reaction.

Real-time quantitative PCR. Real-time quantitative PCR (RT-qPCR) was performed using Fast SYBR Green master mix (Applied Biosystems) and an Applied Biosystems StepOnePlus machine as described previously (42). The synthesized cDNA served as the template. Primer sequences were designed using the National Center for Biotechnology Information Primer-BLAST online tool and were purchased from MilliporeSigma (a full list can be found in Supplemental Table 9). The RT-qPCR was performed in a volume of 20 μL containing 5 μL of cDNA template, 10 μL of 2× Fast SYBR Green master mix (Applied Biosystems), 3 μL water, and 1 μL of the sense and antisense primer (stock concentration of 10 μM each). Mouse or human ribosomal P0 (RPLP0) mRNA served as an endogenous control. The RT-qPCR was performed in duplicate, and the analysis was carried out with the formula 2^-ΔCt, with ΔCt = Ct target gene – Ct endogenous control.

Cell culture of endothelial cells. HUVECs were purchased from Lonza and cultured with endothelial basal medium (EBM, Lonza) supplemented with 10 % FCS (Invitrogen), Amphotericin-B, ascorbic acid, bovine brain extract, endothelial growth factor, gentamycin sulfate, and hydrocortisone (EGM-singleQuots, Lonza) at 37°C and 5% CO₂, with humidified atmosphere. Isolated cardiac mouse fibroblasts were cultured with DMEM high glucose GlutaMAX (Invitrogen) supplemented with 10% FCS (Invitrogen) and penicillin/streptomycin (Roche) at 37°C and 5% CO₂ with a humidified atmosphere. Cell detachment was performed with 0.25% trypsin (Life Technologies, Thermo Fisher Scientific), which was incubated for 2 minutes at 37°C and 5% CO₂ and neutralized with fully supplemented medium. The cell number was determined using the cell counter NucleoCounter (ChemoMetec), as described in the manufacturer’s protocol.

Fibroblast-conditioned medium. To obtain fibroblast-conditioned medium for endothelial cell/fibroblast interaction studies, isolated cardiac fibroblasts from young and old mice were cultured 24 hours after cell isolation in 12-well cell culture plates (Greiner Bio-One GmbH) with 1.5 mL of fully supplemented EBM per well at 37°C and 5% CO₂ with humidified atmosphere. After 24 hours, cell culture supernatants were centrifuged (500 g, 5 minutes) to deplete dead cells and stored at –80°C.

Tube formation assay. In vitro tube formation was performed with 1.5 × 10⁵ HUVECs cultured in 1 mL of fibroblast-conditioned EBM in 12-well plates or with 2 × 10⁴ HUVECs cultured in 100 μL of conditioned medium in 96-well plates (Greiner Bio-One GmbH) that had been coated with 200 or 100 μL Matrigel (356234, Corning). Tube formation was determined after 24 hours by measuring the accumulated tube length in 5 randomly chosen microscopic fields with a computer-assisted microscope using Axiovision 4.5 (Zeiss).

For rescue experiments, medium of aged fibroblasts was supplemented with 4 μg of a serpin-neutralizing (AF3828, R&D Systems) or Efemp1-neutralizing (LS-B11500, LifeSpan Biosciences) antibody.

For the recombinant serpin study, HUVECs were cultured in medium supplemented with 10 ng/mL serpin E1 (CSB-EP021081MO, Cusabio) or serpin E2 (2175-PI-010, R&D Systems) for 24 hours.

Inflammation study. We seeded 80,000 HUVECs in 24-well plates and cultured them in conditioned medium of fibroblasts derived from young and aged mouse hearts in the presence and absence of 10 ng/mL TNF-α. After 24 hours, IL-6 expression was determined using RT-qPCR.

Autophagy. We seeded 40,000 HUVECs in 8-well chamber slides that had been coated with fibronectin before. HUVECs were cultured for 24 hours and stained as follows: HUVECs were fixed for 10 minutes with 4% HistoFix (Carl Roth GmbH + Co. KG), permeabilized with 0.1% Triton X-100, and blocked in 10% donkey serum for 30 minutes. Primary mouse anti-p62 (Novus Biologicals, H00008878-M01) and rabbit anti-LC3 (Novus Biologicals, NB100-2220SS) antibodies were incubated in 10% donkey serum overnight at 4°C. DAPI (1:1000) as well as the secondary anti-rabbit Alexa Fluor 488 (Invitrogen, Thermo Fisher Scientific, A21206) and anti-mouse Alexa Fluor 568 (Invitrogen, Thermo Fisher Scientific, 498389) antibodies were incubated for 1 hour. Autophagosomes were identified as p62/LC3–double positive areas and counted in 4 images per condition.

Proliferation. We seeded 40,000 HUVECs in 8-well chamber slides that were coated with fibronectin. HUVECs were cultured for 24 hours and stained with DAPI. Cell were counted in 4 randomly chosen images.

Secretome analysis. Secreted proteins were detected in supernatants of isolated fibroblasts, which were cultured for 24 hours in DMEM without any supplements. Supernatants were centrifuged (700 g) 5 minutes) to remove cell debris and stored at –80°C. Cytokines and growth factors were detected using the Proteome Profiler Mouse Angiogenesis Array (ARY015, R&D Systems) according to the manufacturer’s protocol.

Picrosirius red staining. Picrosirius red staining was used to determine collagen deposition and fibrosis in paraffin sections of young versus aged mice hearts. In a first step, the paraffin sections were deparaffinized with xylene (twice 10 minutes) and an ethanol series of 100%, 95%, 80%, 70%, and 50% ethanol (5 minutes each step), before washing the sections with water and PBS (5 minutes each step). In a second step, a 0.1% Picrosirius red solution was prepared by dissolving 0.5 g Sirius red (Waldeck GmbH) in 500 mL picric acid (MilliporeSigma) that incubated with the deparaffinized sections for 1 hour. After washing 2 times with acidified water, the sections were dehydrated with 100% ethanol, cleared with xylene, and mounted with Pertex (Medite).

Von Kossa staining. Cardiac calcification was detected on paraffin sections (10 μm) of young (12 weeks) and aged (18–20 months) mouse hearts by using the Silver plating kit acc. to von Kossa (1.00362.0001, Merck). In a first step, the paraffin sections were deparaffinized with xylene (twice 10 minutes) and an ethanol series of 100%, 95%, 80%, 70%, and 50% ethanol (5 minutes each step), before washing the sections with water and PBS (5 minutes each step). In a second step, silver nitrate solution incubated for 20 minutes with the sections. After washing with water, sodium thiosulfate solution was used and incubated for 5 minutes with the sections. Sections were washed with water and nuclei were stained with Kernechtrot (Merck). After washing again with water, sections were dehydrated with an ascending ethanol series (70%, 96%, 100%: 1 minute each), cleared with xylene, and mounted with Pertex (Medite).

Immunostaining of PDGFR-α and osteocalcin. To support single-nucleus RNA-sequencing clustering, fluorescence immunostaining was performed in young (12 weeks) and aged (18–20 months) mouse heart paraffin sections. Hearts were fixed immediately in 4% HistoFix overnight, embedded in paraffin, sectioned at 4-μm thickness using a microtome, and placed on adhesive glass slides. Next, slides were incubated at 60°C for 1 hour and deparaffinized with xylene (twice 10 minutes) and an ethanol series of 100%, 95%, 80%, 70%, and 50% ethanol (5 minutes each step), before washing the sections with water and PBS (5 minutes each step). Sections were unmasked with 0.01 M citrate buffer (pH = 6), washed for 5 minutes with PBS and 0.1% Triton X-100, and blocked for 1 hour with PBS containing 1% BSA, 2% donkey serum, and 0.1% Triton X-100. The primary antibodies anti-osteocalcin (ab93876, Abcam) and anti-Pdgfra (AF1062, R&D Systems) were diluted 1:100 in Pblec buffer (1 mM MgCl₂, 1 mM CaCl₂, 0.1 mM MnCl₂, 1% Triton X-100 in PBS; pH = 6.9) and incubated at 4°C overnight (or for 1 hour at room temperature) on the sections. Sections were washed 3 times with PBS, 0.1% Triton X-100, and the secondary antibodies donkey anti-rabbit Alexa Fluor 488 (A21206, Invitrogen) and donkey anti-goat Alexa Fluor 647 (A21447, Invitrogen) and incubated for 1 hour at room temperature (dilution: 1:300) on the slides. After washing slides again 3 times with 0.1% Triton X-100, slides were mounted with mounting medium containing Hoechst 33342.

Immunostaining for capillary density. Capillary density was determined in 4-μm thick paraffin-embedded sections from young (12 weeks) and aged (18–20 months) mouse hearts. Sections were dried overnight at 37°C, incubated at 60°C for 1 hour, and deparaffinized with xylene (twice 10 minutes). An ethanol series of 100%, 95%, 80%, 70%, and 50% ethanol (5 minutes each step) was applied, before washing the sections with water. Sections were unmasked with 0.01 M citrate buffer (pH = 6), washed for 5 minutes with PBS and 0.1% Triton X-100, and blocked for 1 hour with PBS containing 1% BSA and 0.1% Triton X-100. Biotinylated Isolectin B4 (Vector Laboratories, Maravai Life Sciences, B-1205) was diluted 1:50 and incubated overnight at 4°C. Next, 1:200 streptavidin Alexa Fluor 488 (Invitrogen, Thermo Fisher Scientific, S32354) and 1:400 wheat germ agglutinin Alexa Fluor 647 (Molecular Probes, Thermo Fisher Scientific, W32466) incubated for 1 hour at room temperature. Nuclei were stained with Hoechst 33342. Capillary density was determined as Isolectin B4–positive area per total area in 4 images per condition.

Data availability. The summary of all statistics for each sample can be found in Supplemental Table 1. The full raw and processed data were submitted to the ArrayExpress database of the European Bioinformatics Institute and may be found through accession number E-MTAB-7869.

Statistics. Data of targeted molecular and cell biological procedures described above are represented as means and error bars indicate SEM. After passing Gaussian distribution as assessed by Kolmogorov-Smirnov or D’Agostino-Pearson normality test, statistical power was determined, using 2-sided, unpaired t test for 2-group comparison. For comparisons of more than 2 groups, multiple-group ANOVA with a post hoc Bonferroni’s test was used. (*P < 0.05; **P < 0.01; ***P < 0.001; and ****P < 0.0001.)

Study approval. All animal experiments were approved by the Regional Board of the State of Hessen, Germany.