Single-cell expression profiling reveals dynamic flux of cardiac stromal, vascular and immune cells in health and injury

Abstract

Besides cardiomyocytes (CM), the heart contains numerous interstitial cell types which play key roles in heart repair, regeneration and disease, including fibroblast, vascular and immune cells. However, a comprehensive understanding of this interactive cell community is lacking. We performed single-cell RNA-sequencing of the total non-CM fraction and enriched (Pdgfra-GFP⁺) fibroblast lineage cells from murine hearts at days 3 and 7 post-sham or myocardial infarction (MI) surgery. Clustering of >30,000 single cells identified >30 populations representing nine cell lineages, including a previously undescribed fibroblast lineage trajectory present in both sham and MI hearts leading to a uniquely activated cell state defined in part by a strong anti-WNT transcriptome signature. We also uncovered novel myofibroblast subtypes expressing either pro-fibrotic or anti-fibrotic signatures. Our data highlight non-linear dynamics in myeloid and fibroblast lineages after cardiac injury, and provide an entry point for deeper analysis of cardiac homeostasis, inflammation, fibrosis, repair and regeneration.

Keywords: cell biology; computational biology; heart; mouse; myocardial infarction; scRNA-seq; systems biology.

Figure 1.. TIP scRNA-seq.

(A) t-SNE plots showing detected lineages and sub-populations in TIP across conditions. (B) t-SNE plot of aggregate TIP cells with identified sub-populations. (C) Dendrogram of sub-populations according to average RNA expression. (D) Expression of select marker genes across TIP cells as visualized on t-SNE plots. (E) Cell population percentages across conditions determined to be significantly modulated according to Differential Proporation Analysis (DPA) (p<0.01).

Figure 1—figure supplement 1.. Experimental procedures, population proportions and gene expression characterstics of sub-populations within TIP scRNA-seq.

(A) Diagram of experimental procedures. (B) Percentage of total cells in each cell-population according to experimental condition. (C) Percentage of total cells in clusters as grouped together according to high-level cell type and as according to experimental condition. (D) Dot plot showing expression of top upregulated genes across TIP populations. (E, F) Box plot showing number of UMIs (E) and genes (F) detected across TIP sub-populations after quality control filtering.

Figure 1—figure supplement 2.. Gene expression data for TIP sub-populations.

(A) Distribution of Pdgfra and Pdgfra-GFP expression across cell populations in TIP. (B) Expression of selected endothelial markers as visualized in box and t-SNE plots. (C) Expression of selected cell cycle marker genes on t-SNE plots. (D) Dot plot visualization according to condition of top marker genes for lymphoid lineage populations.

Figure 1—figure supplement 3.. Differential proportion analysis (DPA) procedure and evaluation.

(A,B) Procedure for performing DPA. (C) Percentage of cell populations between sham conditions from GFP⁺ cohort. (D) Comparisons of p-values from Fisher’s exact test and DPA for evaluating proportion changes between sham experiments. (E) Comparisons of specificity for Fisher’s exact test and DPA on control replicate simulation experiment. (F–H) Sensitivity, specificity and precision for simulated control vs condition experiments with comparisons between Fisher’s exact test and DPA.

Figure 1—figure supplement 4.. Clustering of TIP scRNA-seq prior to removal of minor hybrid populations.

(A) t-SNE plot of aggregate TIP cells with all identified 29 sub-populations prior to filtering of minor hybrid populations. (B) Dendrogram of all 29 sub-populations according to average RNA expression. (C) Expression of example marker genes that illustrate hybrid populations (F-EC, M2MΦ-EC, EC-L1, EC-L2 and BC-TC) as visualized in box and t-SNE plots. (E) Percentage of total cells in each sub-population according to experimental condition.

Figure 1—figure supplement 5.. A typical workflow of sequential gating strategy for doublet exclusion is shown.

(A) Scatter gate, FSC-A vs SSC-A, was plotted to exclude cell debris; the three subsequent ‘doublet discrimination’ gating strategies to exclude doublets are: (B) FSC-H vs FSC-A plot, (C) FSC-H vs FSC-W and (D) SSC-H vs SSC-W. Blue and red gates depict ‘single cell’ and ‘doublet’ events, respectively. Only the cells in the blue gates were sorted and/or analyzed. The scatter profiles of putative doublets (red dots) are also illustrated in panel (C) and (D).

Figure 1—figure supplement 6.. FACS data supporting F-EC population.

FACS plots showing expression of GFP or CD31 (axis label) in the (A) wild type sample stained with isotype control antibody; (B) wild type sample stained with anti-CD31 antibody; (C) Pdfgra^GFP/+ sample stained with isotype control antibody; (D, E) Pdfgra^GFP/+ sample stained with anti-CD31 antibody 7 days after Sham (D) or MI (E) surgery. Numbers indicate percentage of GFP⁺CD31⁺ double positive cells in live single cells in respective plot (mean numbers indicated).

Figure 1—figure supplement 7.. FACS data supporting the M2 MΦ -EC population.

FACS plots showing expression of indicated antigen (axis label) in the (A) wild type sample stained with isotype control antibodies and Pdfgra^GFP/+ sample stained with indicated antibodies 7 days after Sham (B) or MI (C) surgery. Numbers indicate percentage of gated cells in live single cells in respective analyses (mean numbers indicated).

Figure 1—figure supplement 8.. Atrial Fibrilation (AF) associated genes in TIP.

Overlap of genes associated with AF from Roselli et al. (2018) and genes upregulated in TIP sub-populations with log2 fold-change >1. Displayed are genes upregulated in EC (top) or mural (bottom) cell sub-populations with expression visualized in box and t-SNE plots. Key indicates the location of sub-populations on the t-SNE plot.

Figure 2.. Cardiac Mo/MΦ populations.

(A) t-SNE plot showing extracted Mo/MΦ populations. (B) Expression of select immune cell markers as visualized in box-plots and t-SNE plots. Arrows indicate Ccr2^low sub-population of M2MΦ. (C) Dot-plot of top upregulated genes for each Mo/MΦ population where color indicates experimental conditions. (D) Heatmap of differentially expressed genes between Mo/MΦ populations with representative significant GO Biological Process terms.

Figure 2—figure supplement 1.. Mo and MΦ marker genes and Diffusion Map analysis.

(A) Expression of selected genes marking Mo and MΦ sub-populations visualized using box and t-SNE plots. (B) Expression of H2-Aa (MHC-II) and Igf1 between sub-populations of M2MΦ defined according to Ccr2 expression (Ccr2^high and Ccr2^low) as visualized in box-plots. (C) Major Mo/MΦ populations visualized on diffusion components 1–3.

Figure 3.. Cell-cell ligand-receptor network analysis.

(A) Hierarchical network diagram of significant cell-cell interaction pathways. Arrows and edge color indicates direction (ligand:receptor) and edge thickness indicates the sum of weighted paths between populations. (B) Comparison of total incoming path weights vs total outgoing path weights across populations. (C) Summed ligand weights across souce ligand and receptor target paths for top ligands in MYO. (D) Tree plot showing outgoing connections from the Glial cells. Top node refers to source population, second layer to ligands, third layer to receptors and leaf nodes represent target populations.

Figure 3—figure supplement 1.. Pdgfra-GFP⁺ cells localization in healthy and diseased hearts.

(A–B) Representative image of GFP (green) and CD31 (red) co-immunostainings in sham (A) and MI-day 3 (B) hearts. Arrowheads show GFP⁺ cells in close proximity/contact with CD31⁺ cells. Scale bar - 20 μm.

Figure 4.. Pdgfa-GFP⁺ scRNA-seq.

(A) t-SNE plot of GFP⁺ cells separated according to experimental condition (sham, MI-day 3, MI-day 7). (B) t-SNE plot showing aggregate of GFP⁺ cells across conditions. (C) Dendrogram of populations determined by average RNA expression in populations. (D) Percentages of cells in each population according to experimental condition. Stars indicate significant change across conditions according to DPA (p<0.01). (E) Expression of select genes in different populations as visualized in box and t-SNE plots. (F) Dot-plot of top five upregulated genes for each population where color indicates strength of expression and size of dot represents percentage of cells expressing the gene.

Figure 4—figure supplement 1.. FACS profiles and scRNA-seq analysis for GFP⁺/CD31^- cells.

(A) FACS plots and gating strategies for sorting GFP⁺/CD31^- population at 7 days post-sham or MI surgery. (B) FACS plots showing SCA1 expression profiles in GFP⁺/CD31^- cells at 7 days post-sham or MI surgery. (C) Boxplot showing number of UMIs detected per-cell across experimental conditions in GFP⁺/CD31^- fraction after quality control filtering. (D) Boxplot showing number of genes detected per-cell across experimental conditions in GFP⁺/CD31^- fraction after quality control filtering. (E) t-SNE visualization of scRNA-seq clusters for all sham and MI time-points in GFP⁺/CD31^- cells. (F) Cell population proportions across all experimental conditions.

Figure 4—figure supplement 2.. Comparison of GFP⁺ populations with FACS-sorted S⁺P⁺ cells.

(A) ROC curves showing prediction accuracy of iRF classifier for predicting populations in sham cells. ROC calculated across a 10-fold cross-validation test. (B) Counts of iRF predicted cell identities across two scRNA-seq (Fluidigm) experiments on S⁺P⁺ cells. (C) Thy1 expression across conditions. (D) Distribution of Thy1.2 (CD90.2) protein expression by FACS in S⁺Pdgfra-GFP⁺cells. (E) Colony counts comparing S⁺Pdgfra-GFP⁺ with Thy1.2 fractions (low vs high) isolated by FACS (n = 3). Indicated p-values are derived from Student’s t-test.

Figure 4—figure supplement 3.. Expression of selected marker genes across GFP⁺ populations visualized in box and t-SNE plots.

Circles indicate relevant population for the displayed marker.

Figure 4—figure supplement 4.. Population-specific expression of transcription factors marking fibroblast and myofibroblast sub-populations.

Figure 4—figure supplement 5.. Comparisons with Skelly et al. (2018) and Gladka et al. (2018) scRNA-seq data-sets.

(A–C) Reanalysis of Skelly et al. scRNA-seq data. (A) t-SNE plot showing main fibroblast populations detected from clustering analysis and identities according to those described in GFP⁺. (B) iRF predictions of GFP⁺ sham populations on fibroblasts from Skelly et al. data-set. (C) Expression of representative marker genes according to box plots and as visualized on t-SNE plots. Shows fibroblast markers (Pdgfra, Ddr2, Col1a1) and sub-population-specific markers. (D) Expression of Ckap4 gene, which was identified by Gladka et al. as a novel marker of activated fibroblasts, as visualized in box and t-SNE plots across TIP and GFP⁺ sub populations.

Figure 5.. Features of the F-WntX population.

(A) Differentially expressed genes in F-WntX overlaid on Wnt pathway maps. Wnt node includes genes Wnt5a and Wnt16. (B) Detection rate, representing the percentage of cells expressing a gene, across all conditions for cells in F-WntX or all GFP⁺ cells combined. (C) Example GO BP terms over-represented (FDR < 0.05) in genes upregulated in F-WntX compared to F-SL/F-SH populations. (D) Tree plot showing ligand-receptor connections from F-WntX to EC sub-populations as calculated in TIP. Top node refers to source population, second layer to ligands, third layer to receptors and leaf nodes represent target populations. (E) Examples of F-WntX:EC ligand and corresponding receptor expression as visualised in box and t-SNE plots. For each ligand the corresponding receptor is immediately below.

Figure 5—figure supplement 1.. Activation and paracrine signature of F-WntX cells.

(A) Expression of markers of activation visualized on condition-specific t-SNE plots. (B) GO Molecular Function terms over-represented (FDR < 0.05) in genes upregulated in F-WntX compared to F-SL/F-SH. (C, D) Venn diagram of overlap between F-WntX marker genes in ‘Regulation of cell proliferation’ category and MF categories ‘Signaling receptor binding’ (C) and ‘Growth factor binding’ (D). (E) Examples of ligands upregulated in F-WntX and corresponding receptors as expressed in TIP and visualized in box and t-SNE plots. For each ligand the corresponding receptor is immediately below. Circle indicates the location of F-WntX on the t-SNE plot.

Figure 6.. WIF1 localization and co-expression in injured and uninjured hearts of Pdgfra-GFP⁺ mice.

(A) Representative image of WIF1 (red), GFP (green) and Wheat Germ Agglutinin (WGA, grey) co-immunostaining showing the border zone at MI-day 3. Arrowheads show WIF1⁺ cells. Scale bar - 50 μm. (B–C) Representative images of WIF1 (red) and GFP (green) co-immunostainings showing left ventricle (sham, (B) or infarcted border zone at MI-day 3 (C). Arrowheads show WIF1⁺ cells, Asterix shows WIF1⁺GFP⁺ cells. Scale bars - 20 μm. (D–G) Representative images of co-immunostainings for WIF1 (red), GFP (green) and markers (gray) for golgi (GM130, D, Arrowheads show WIF1⁺GM130⁺ cells), proliferation (Ki67, E, Arrowheads show WIF1⁺Ki67⁺ cells), smooth muscle cells and myofibroblasts (α-SMA, F, (F’) showing an example of a WIF1⁺GFP⁺α-SMA⁺ cell from another section), and leukocytes (CD45, G, Arrowheads show WIF1⁺CD45⁺ cells, Asterix shows GFP⁺CD45⁺ cells). Scale bar - 20 μm. (H) Representative image of co-immunostaining for WIF1 (green), WGA (gray) and endothelial cell marker CD31 (green). Arrowheads show WIF1⁺ cells in close proximity/contact with CD31⁺ cells. Scale bar - 20 μm. (I) Quantification of marker-positive cells in the infarcted border zone of MI-day 3 hearts. n = 4.

Figure 6—figure supplement 1.. WIF1 expression pattern in E14.5 embryos.

Representatives images of sagittal sections of E14.5 embryos immunostained with WIF1 (red) antibody. Scale bar is 50 μm.

Figure 6—figure supplement 2.. WIF1 protein expression after sham or MI.

(A–D) Representative images of WIF1 (red) and GFP (green) co-immunostainings in the left ventricle (sham, A) or in the infarcted border zone of hearts, 1 (B), 3 (C) or 7 (D) days after surgery. Arrowheads show WIF1⁺ cells, Asterix shows WIF1⁺GFP⁺ cell. Scale bar - 20 μm.

Figure 7.. Diffusion Map analysis of GFP⁺ cells.

(A) 3D Diffusion Map of main fibroblast/myofibroblast populations with cells colored according to population. (B) 2D Diffusion Map facetted according to experimental condition. (C) Expression of marker genes on main trajectories of diffusion components across conditions. (D) Heatmap of differentially expressed genes with representative GO Biological Process terms.

Figure 8.. Time-point-specific analysis of GFP⁺ scRNA-seq.

(A,B) t-SNE visualization of GFP⁺ populations 3 days post sham/MI (A) and 7 days post sham/MI (B). (C,D) 3D Diffusion Map analysis of day 3 major populations (C) and day 7 major populations (D). (E) Heatmap of upregulated genes in day 3 injury-response populations. (F) Heatmap of differentially expressed genes between myofibroblast sub-populations. (G) Gene expression visualized in box and t-SNE plots for myofibroblast sub-population marker genes.

Figure 8—figure supplement 1.. Reanalysis of GFP-day 3 and GFP-day 7 data-sets separately.

(A) Expression of F-CI marker genes on GFP-day 3 dot plot (B) AUC, sensitivity, specificity and precision for iRF classifier trained to predict GFP-day 3 cell populations after 10-fold cross-validation. (C) iRF prediction scores for GFP-day 7 cells as visualized on t-SNE plots. (D) Counts of GFP-day 7 cells predicted to correspond to GFP-day 3 populations as determined by iRF prediction (score >0.5). (E) Percentage of GFP-day 3 F-Cyc cells reassigned to alternative populations after removing cell cycle genes. (F) Percentage of GFP-day 7 F-Cyc cells reassigned to alternative populations after removing cell cycle genes. (G) 2D Diffusion Map plot of F-Act, MYO-1, MYO-2 and MYO-3 from the GFP-day 7 data-set.

Figure 9.. Schematic summary of the flux and pseudotime differentiation dynamics of GFP⁺ populations between sham, MI-day 3 and MI-day 7.

Populations are ordered in pseudotime from unactivated (top) to most activated/mature (bottom). Arrows connecting populations indicate direction of proposed differentiation/pseudotime trajectory. Colored arrows indicate whether the population appears to expand or diminish relative to the previous time-point.