Single

Clinical Perspective

What Is New?

• Our study identifies a new subpopulation of cardiac fibroblasts (CFs), reparative CFs, characterized by a distinct transcriptional profile, including Cthrc1, with a major role in the fibrotic healing response after myocardial infarction (MI).

• Experiments performed in a mouse model deficient in Cthrc1 revealed an increase in mortality after MI and a decrease in the fibrotic response.

• We identified CTHRC1⁺ CFs in a swine model of MI and in patients with MI and dilated cardiomyopathy, supporting the conservation of the reparative CF subpopulation in humans.

What Are the Clinical Implications?

• Understanding the role of CFs after MI is highly relevant to identify new therapeutic targets for controlling adverse remodeling after MI.

• The present study provides a unique resource for the cardiovascular field by characterizing the subpopulation of CFs responsible for cardiac repair after cardiac fibrosis and identifying molecular mechanisms implicated in the generation of this subpopulation.

• The conservation of this population in swine and humans and the relation between Cthrc1 expression and cardiac function after MI suggest that this molecule might be a potential target in patients with MI.

Introduction

Cardiac fibroblasts (CFs) represent only 10% of the total number of cells in the myocardium; however, they play a critical role in the structural and mechanical maintenance of the heart.^1,2 After myocardial infarction (MI), CFs become activated, orchestrating a fibrotic response that leads to the generation of a collagen scar, which prevents cardiac rupture.^3,4 The absence of reliable markers has traditionally hampered our understanding of the specific role of CFs in cardiac homeostasis. Several studies using lineage-tracing reporter strains have allowed us to track the origin and roles of fibroblasts, showing that the fibroblast response to cardiac injury is heterogeneous (reviewed in reference 3). These studies revealed that different CF subtypes may play different roles during the healing process that follows MI.^5–7 In this context, Fu et al⁸ have recently described a new subpopulation of activated CFs, the matrifibrocytes, which support the mature scar and are characterized by the expression of extracellular matrix (ECM) and tendon genes.

This highlights the need for a better understanding of CF heterogeneity and its effect on processes that mediate repair of the ischemic myocardium. The development of single-cell RNA sequencing (scRNA-seq) represents an ideal tool to address this knowledge gap. Two recent publications have revealed the presence of 3 subpopulations of CFs in pathologic conditions at single-cell resolution using a sorting and robot-assisted transcriptome sequencing protocol.^7,9 In accordance with previous studies, both articles show that activated CFs are characterized by the expression of Postn. In contrast, Kretzschmar et al¹⁰ described a total of 11 subpopulations of CF after pooling samples from neonatal, adult, and pathologic hearts. All these studies have poor representation of CFs because of the small number of collected cells (63 CFs/185 cardiac cells in reference ⁷, 935 cardiac cells in reference ⁹, and 243 CFs in reference ¹⁰). In contrast, Skelly et al¹¹ performed an unbiased analysis of 10 519 cardiac cells distinguishing 4 subpopulations of CF but only during homeostasis. More recently, Farbehi et al¹² described 8 CF subpopulations within a total of 11 that emerged from unbiased clustering of 16 787 Pdgfrα-positive cardiac interstitial cells isolated from healthy hearts and hearts 3 and 7 days after infarction. No scRNA-seq study to date has focused exclusively on CFs, their anatomic location, and their dynamics during the course of MI.

We used scRNA-seq to identify a new subpopulation of CFs that mediates ventricular remodeling associated with different cardiac fibrotic processes. This subpopulation is characterized by a specific transcriptomic signature that includes Cthrc1 (collagen triple helix repeat containing 1), an essential molecule involved in the synthesis and deposition of the fibrotic scar.^13,14 We describe and characterize the role of some putative regulatory transcription factors (TFs) and validate the noncanonical transforming growth factor (TGF)–β signaling pathway as a mediator of this specific signature. Using a knockout model for genetic ablation of CTHRC1, we provide evidence for the essential role of this specific fibroblast population during cardiac healing. We provide evidence of the conservation of this CTHRC1⁺ CF subpopulation in a preclinical swine model of myocardial ischemia and in patients with different types of cardiac fibrosis.

Methods

All data, analytic methods, and study materials are available from the corresponding author on request.

Animal Models

All animal requisitions, housing, treatments, and procedures were performed according to all state and institutional laws, guidelines, and regulations. All studies were approved by the Ethics Committee for Animal Research at the University of Navarra and the Government of Navarra. Studies conducted at Maine Medical Center Research Institute were approved by its Institutional Animal Care and Use Committee and were compliant with The Guide for the Care and Use of Laboratory Animals.

Human Samples

The study protocol was approved by the Medical Ethics Committee and informed consent was obtained from all patients. For RNA studies, human heart samples were provided by Professor S. Janssens. For histologic analysis, human heart samples were provided by the Maine Medical Center BioBank.

Statistical Analysis

In vitro scratch and morphometric assay experiment statistical significance was analyzed using a nonparametric 1-way analysis of variance with a Kruskal-Wallis post hoc test. For the Kaplan-Meier survival curve, log-rank Mantel-Cox test was used to determine statistical difference between the survival curves of the 2 groups of animals. Quantification of collagen deposition statistical significance was analyzed by an unpaired t test.

Statistical analysis to obtain the significance when comparing Z score distributions was done with a 2-sided Wilcoxon rank sum test. Differential expression analysis for normalized values was performed using negative binomial generalized linear models.

A 2-way analysis of variance was used to determine statistical significance between the relative quantitative polymerase chain reaction values in different anatomic regions for POSTN, COL1α1, and CTHRC1 in pig samples. For comparison of the echocardiography, resonance, and infarct size between reduced and preserved pigs, a 2-tailed t test with a Mann-Whitney post hoc test was used.

Data Resources

All scRNA-seq, bulk RNA sequencing, and ATAC-seq (assay for transposase-accessible chromatin sequencing) data are available in a SuperSeries at the National Center for Biotechnology Information Sequence Read Archive database under accession number GSE132146.

Detailed methods are provided in the Data Supplement, including Supplemental Methods, Tables I–VI, Figures I–XVI, Movie I, and File.

Results

Characterization of CF Heterogeneity

To follow the dynamics of the global CF population during MI, we used Col1α1-GFP reporter mice. This strain has been shown to homogeneously label fibroblasts from different origins¹⁵ during homeostasis and in different cardiac fibrosis models.^1,6,7,16 Histologic analysis of healthy and infarcted myocardium at different time points showed changes in the amount and distribution pattern of GFP (green fluorescent protein)–positive cells around the injured site (Figure 1A). To validate this model, we isolated putative CFs (GFP⁺/CD31⁻/CD45⁻), endothelial cells (GFP⁻/CD31⁺/CD45⁻), and bone marrow–derived cells (GFP⁻/CD31⁻/CD45⁺) from healthy myocardium along with myocardium 7, 14, and 30 days postinfarction (dpi; Figure 1B). No GFP⁺ cells were found in the CD31⁺ or CD45⁺ cell fractions at any time point by fluorescence-activated cell sorting or immunohistochemistry (Figure 1B and Figure IA in the Data Supplement). No transcriptional markers for CFs were found in these populations of cells. The transcriptional profiles of putative CFs and tail/dermal fibroblasts (GFP⁺/CD31⁻/CD45⁻) were shown to be closely similar (Figure IB–ID in the Data Supplement). An increased percentage and total number of GFP⁺ cells was observed in the infarct zone (IZ) and border zone (BZ) after MI (Figure IIB in the Data Supplement). Col1α1-GFP⁺/CD31⁻/CD45⁻ cells expressed a set of surface membrane markers and a transcriptomic profile consistent with classic CFs (Figure IIA and IIC in the Data Supplement) and mEFSK4⁺/CD31⁻/CD45⁻ CFs^1,17 (Figure IID and IIE in the Data Supplement). Taken together, these results indicate that cardiac GFP⁺ cells are bona fide fibroblasts, validating the Col1α1-GFP model to study CF biology during cardiac repair.

• Download figure

• Download PowerPoint

Figure 1. Cardiac fibroblasts (CFs) are a heterogeneous population of cells. A, Overview of healthy hearts and infarcted hearts at 7, 14, and 30 days postinfarction (dpi) using bright field and fluorescence microscopy; overview and detail of transverse sections show GFP⁺ (green fluorescent protein) cell distribution in the free wall of the left ventricle (LV) in healthy hearts and infarcted hearts 7, 14, and 30 dpi in the infarct and remote zones and detail of GFP⁺ cardiac interstitial cells. Cardiac troponin I⁺ (cTN) cardiomyocytes, gray; CD31⁺ endothelial cells, red; nuclei, blue (DAPI [4′,6-diamidino-2-phenylindole]). B, Representative gating for isolation of GFP⁺/CD31⁻/CD45⁻ cardiac interstitial cells as performed in healthy hearts and infarcted hearts at 7, 14, and 30 dpi. From left to right: gating for singlets, metabolically active cells (Calcein⁺/TO-PRO-3⁻), viable cells (GFP⁺/TO-PRO-3⁻), and GFP⁺ cells that are negative for CD45 and CD31. C, Schematic representation of the experimental design. Col1α1-GFP mice were subjected to myocardial infarction and GFP⁺ CFs were isolated at different time points (healthy hearts and infarcted hearts at 7, 14, and 30 dpi). These cells were used for single-cell RNA sequencing (scRNA-seq). The functional role of CFs was analyzed with a mouse knock-out model and their regulatory mechanisms with assay for transposase-accessible chromatin sequencing (ATAC-seq) and in vitro assays to validate the hypothesis generated. D, t-Distributed stochastic neighbor embedding (t-SNE) plot of the global CF population (29 176 cells) comprising the 4 different times. Plot is color-coded by clusters (A through K) identified through unsupervised analysis. Dashed lines delimit clusters. E, Heatmap shows log-normalized unique molecular identifiers (UMIs) for scRNA-seq analyzed cells. Top and side bars indicate clusters. Traditional marker genes for CFs are indicated (left). Dot plot of expression and specificity of top markers (2) for identified clusters. Dot size represents percentage of cells per cluster expressing the given marker and color represents the relative expression. Scale bars, 100 µm (37 µm in top right picture in A). Endo indicates endocardium; Epi, epicardium; and LAD, left anterior descending.

To characterize the fibroblast response after MI at single-cell resolution, we profiled the single-cell transcriptomes of GFP⁺ CFs (Figure 1C). Transcriptomes from each time point (7079 cells in healthy myocardium, 10 448 cells at 7 dpi, 8337 cells at 14 dpi, and 6805 cells at 30 dpi) were subjected to quality control filtering and merged into a single data set of 29 176 individual cells using a canonical correlation approach (Figure III in the Data Supplement). Unsupervised clustering of the global dataset revealed 11 clusters of GFP⁺ cells (Figure 1D). Ten of these clusters (A through J) represent different fibroblast subpopulations, characterized by high levels of fibroblast-associated molecules, whereas cluster K comprises a population with high expression of classic pericyte markers¹¹ (Figure 1E and Figure IVA and IVB in the Data Supplement). Clusters B, D, I, and J showed specific expression profiles, representing potential functional CF subgroups, whereas the remaining clusters (A, C, E, F, G, and H) showed less specific transcriptomic identities, likely reflecting intermediate CF types within a larger population (Figure 1E and Figure IVC in the Data Supplement). We also compared our sorted GFP⁺ CFs with CFs identified in a recent study at scRNA-seq resolution (Figure V in the Data Supplement) and found that all GFP⁺ CFs fall within previously identified CFs.¹²

Cluster B Represents a Subpopulation of Activated Periostin⁺ Fibroblasts Localized in the Infarcted Myocardium

To explore the potential role of the different clusters during MI, we analyzed changes in the percentage of cells in each cluster. Clusters F and H through K remained constant at the different time points (Figure VIA in the Data Supplement), whereas cluster B manifested significant changes after MI. Only 2.3% of the total GFP⁺ cells belonged to cluster B in healthy myocardium and the percentage rose sharply after MI (12% at 7 dpi and 34% at 14 dpi) followed by a decrease at 30 dpi (12%; Figure 2A and 2B and Figure VIA in the Data Supplement). Cluster B’s transcriptomic signature was significantly enriched in pathways and gene ontology terms related to ECM organization, cell proliferation, and cell–substrate adhesion, in comparison with other clusters (Figure 2C and Figure IVD in the Data Supplement). This was reflected by high and specific expression of structural molecules, including Fmod and Comp, and by enzymes involved in collagen metabolism including Ddah1, Lox, and Ptn.^18,19Cthrc1, a gene linked to vascular remodeling and fibrotic processes,^20–24 was identified as a top marker of cluster B (Figure 2D and Figure VIB in the Data Supplement). Changes in the prorepair expression signature of cluster B fibroblasts were observed at 7, 14, and 30 dpi and were consistent with potential functions involved in specific phases of the repair process (Figure 2E and Figure VIC in the Data Supplement). The highly specific expression of ECM-related genes in cluster B suggests a profound association of this population with the formation of the fibrotic scar.

• Download figure

• Download PowerPoint

Figure 2. Dynamics of cardiac fibroblast (CF) heterogeneity reveal a unique subpopulation that responds to myocardial infarction. A, t-Distributed stochastic neighbor embedding (t-SNE) plot shows CF heterogeneity at different time points. Number indicates percentage of cells in cluster B. B, Representation shows the dynamics of the percentages of cells by cluster along infarction. Line color-coded by clusters (A through K). C, Network representation of enriched pathways based on cluster B markers. Dot sizes represent number of cluster B markers annotated for each pathway and color scale statistical significance for each function. D, t-SNE representation of normalized expression of Cthrc1, Ddah1, Lox, Comp, Fmod, and Ptn comprising the 4 time points. Dashed lines delimit cluster B. E, Scaled gene expression heatmap shows transcriptional dynamics of cluster B along time points. F, Representative immunohistochemistries of Col1α1-GFP infarcted hearts at 14 days postinfarction (dpi) show the spatial location of CTHRC1, DDAH1, or FMOD (red) in the infarct zone. GFP⁺ (green fluorescent protein), green; nuclei (DAPI [4′,6-diamidino-2-phenylindole]), blue; colocalizations, yellow (arrows). Scale bars, 100 µm. G, Experimental design: CFs (GFP⁺/CD31⁻/CD45⁻) from remote, border, and infarct zones were sorted for bulk RNA sequencing (top). Boxplot shows the Z score distribution of cluster markers in the 3 zones at 7 dpi (Wilcoxon signed rank test, n=2; middle). Bar plot shows normalized expression values for top cluster B markers (bottom; linear model differential expression, n=2). Data are mean±SD. Scale bars, 100 µm. *P≤0.05, **P≤0.01, ***P≤0.001. ADAM indicates A disintegrin and metalloproteinase; and TGF, transforming growth factor.

Immunohistochemical analysis of 3 specific cluster B markers (Cthrc1, Ddah1, and Fmod) revealed that these cells were almost exclusively located in the IZ and BZ at 7, 14, and 30 dpi (Figure 2F and Figure VII in the Data Supplement). In agreement, zonal RNA sequencing profiling of the infarcted heart at 7 dpi showed a clear enrichment in the cluster B signature in the infarct area compared with the remaining clusters (Figure 2G). Taken together, our results identify cluster B as a defined subpopulation of CFs that emerges in response to MI, presents a strong and dynamic profibrotic signature, and localizes to the damaged tissue.

To analyze the origin of cluster B, both RNA velocity and latent time analysis were performed. Using RNA velocity,²⁵ we observed a dynamic movement among all clusters of CFs in healthy heart (Figure VIII in the Data Supplement). At 7 dpi, the previously identified transitions were clearly reduced, and a specific transition from cluster F to cluster B was observed. This transition was prominent at 14 dpi, suggesting that cluster B cells were likely originating from cluster F. At 30 dpi, most of the identified transitions were reduced, including F to B, resembling the dynamics between subpopulations observed in healthy myocardium (Figure VIIIA in the Data Supplement). Using latent time analysis,²⁶ cluster B cells grouped in regions with the highest latent time values. The expression of Cthrc1 was almost exclusive to cluster B, whereas other marker genes for quiescent (Cd90/Thy1 or Pdgfrα/Cd140a) or activated CFs (periostin [Postn] or Acta2/αSMA) were expressed in earlier stages of activation (Figure VIIIB in the Data Supplement). This pattern suggests that cluster B fibroblasts correspond to the final activation stage of a subset of Postn⁺ CFs generated in response to MI. In this process, expression of Cthrc1 seems to be specific for the final stages of activation and associated with the healing process through scar formation.

To determine the differences between cluster B cells and other subpopulations of activated CFs,^4,8 we compared the transcriptomic profile of CFs that express or coexpress the specific markers Postn (activated CFs), Acta2 (matrifibrocites), or Cthrc1 (cluster B fibroblasts). Depending on the combined expression of these 3 markers, different subsets of activated CFs appeared (Figure IXA in the Data Supplement). Again, Cthrc1⁺ CFs emerged as a subset of Postn⁺ CFs, but with enrichment in gene ontology terms related to ECM assembly and organization, as well as collagen fibril organization (Figure IXB in the Data Supplement). In agreement with these results and our scRNA-seq analysis, different subpopulations of activated CFs were detected at 7 dpi by immunohistochemistry (Figure IXC in the Data Supplement). A small number of CTHRC1⁺ CFs were identified at 60 dpi (Figure IXD in the Data Supplement).

To determine whether cluster B can be identified in other models of cardiac fibrosis, we examined the presence of CTHRC1⁺ CFs in a model of chronic cardiac fibrosis generated by infusion of angiotensin II for 28 days.⁵ As shown in Figure X in the Data Supplement, CTHRC1⁺ CFs were detected 2 days after the beginning of treatment with angiotensin II (dpt), with an increase at 7 dpt. At 14 dpt, this population decreased and was no longer detected 28 days after infusion of angiotensin II.

Collectively, these results suggest that cluster B cells represent a specific subpopulation of activated Postn⁺ CFs with a major role in reparative cardiac fibrotic processes. Hence, we termed these cells reparative cardiac fibroblasts (RCFs).

Characterization of the Molecular Regulation of the RCF Response

To identify the molecular regulation underpinning the generation of RCFs, we undertook different approaches. First, we leveraged publicly available ChIP-seq (chromatin immunoprecipitation combined with sequencing) datasets to identify TFs whose binding patterns are enriched in the vicinity of RCFs genes. This approach identified several TFs, such as SOX9 and SMAD3 (Figure 3A). Some of the top marker genes of RCFs showed binding motifs for SOX9 (Figure 3B). Sox9 overexpression in cultured CFs induced 23% of the RCFs signatures (28 genes; fold change >1.5, P<0.05; Figure 3C), similar to the response observed after incubation of CFs with TGF-β (33 genes; fold change >1.5, P<0.05), a classic regulator of CF activation^27,28 (Figure 3D). These results suggested the potential of SOX9 to promote an RCF phenotype and are consistent with previous studies describing CF activation programs.^28,29 However, Sox9 overexpression only partially explained the transcriptional signature of RCFs.

• Download figure

• Download PowerPoint

Figure 3. CD200⁺/CD146⁻ cardiac fibroblasts (CFs) provide the most specific characterization of reparative CFs (RCFs). A, Transcription factor (TF) target gene enrichments. Dot size represents enrichment P value of TF and color represents log 2 transformed expression in CD200⁺ bulk RNA sequencing. B, SOX9 DNA binding motif sequence logo and its location in some of the top markers for cluster B. Small red rectangles below the gene sequence show confirmed SOX9 binding sites, as determined by JASPAR. C and D, Volcano plots show differential gene expression of in vitro grown CFs overexpressing Sox9 (left) or treated with transforming growth factor (TGF)–β (right). Genes with log fold change of ±1.5 and P<0.05 were considered differentially expressed. Red dots represent RCF markers. E, Violin plots show single-cell normalized expression of selected surface markers in the pooled CF population (healthy hearts and infarcted hearts at 7, 14, and 30 days postinfarction [dpi]). F, Transverse sections of Col1α1-GFP hearts at 7 dpi. Scale bars, 1 mm. Immunofluorescence analysis of GFP⁺ (green fluorescent protein; green), CD200⁺ (red), COL1α1 (gray), and DAPI (4′,6-diamidino-2-phenylindole)/nuclei (blue) in healthy left ventricle (LV), and infarct zone (IZ; left) and remote zone (right) at 7 dpi. Colocalization of GFP⁺ and CD200⁺ in yellow (arrows) in IZ (and colocalized with COL1α1 in light yellow). Arrowheads indicate GFP⁺/CD200⁻ cells and asterisks indicate GFP⁻/CD200⁺. Scale bars, 100 µm. Quantification of GFP⁺/CD200⁺ cells in healthy and at 7 dpi (right, top). G, Gating strategy for isolation of GFP⁺/CD31⁻/CD45⁻/CD200⁺/CD146⁻ (CD200⁺) and GFP⁺/CD31⁻/CD45⁻/CD200⁻/CD146⁻ (CD200⁻) cardiac interstitial cells at 7 dpi. H, Boxplot representation of Z score distributions for cluster markers in CD200⁺ and CD200⁻ CF (top), Wilcoxon signed rank test, n=2. Normalized expression bar plots (mean±SD) for top RCF markers in both subpopulations (bottom), linear model differential expression, n=2. Data are mean±SD. *P≤0.05, **P≤0.01, ***P≤0.001. FACS indicates fluorescence-activated cell sorting.

To further identify TFs involved in the generation of RCFs, we aimed to profile chromatin accessibility (ATAC-seq) at 7 dpi. We examined the expression patterns of membrane surface markers in our single-cell dataset to identify a specific combination that would allow the isolation of RCFs. None of the classic markers of fibroblasts such as Thy1 (Cd90), Pdgfrα (Cd140a), or Pgp1 (Cd44)^17,30 was specific for RCFs (Figure 3E). Only Ox-2/Cd200 was expressed in 62.8% of cluster B cells and in 50% of cluster K cells, which comprise RCFs and pericytes, respectively. The percentages of Cd200 were substantially lower in other clusters (Figure XIA in the Data Supplement). CD200⁺/GFP⁺ cells colocalized with COL1α1 protein in the center of the scar (Figure 3F and Figure XIC in the Data Supplement). We next profiled fluorescence-activated cell sorted CD200⁻ and CD200⁺ fibroblasts at 7 dpi after negative selection for CD146 (Mcam) to discard pericytes (Figure 3G). CD200⁺/CD146⁻ CFs showed a significant similarity with the aggregated expression signature of cluster B cells (Figure XIB in the Data Supplement) and significantly higher levels of cluster B, F, G, and H markers in comparison with CD200⁻ CFs (Figure 3H). Using the single-cell data from 7 dpi, 37% of Cd200⁺ cells were assigned to RCFs, whereas other clusters showed smaller proportions (Figure XID in the Data Supplement). No combination of surface markers allowed for the unequivocal selection of cluster B CFs, but the expression of Cd200⁺/Cd146⁻ provided the most specific characterization of RCFs. Based on these results, we performed the chromatin accessibility patterns (ATAC-seq) on CD200⁺/CD146⁻ and CD200⁻/CD146⁻ at 7 dpi, as well as on the GFP⁺ CF population at different time points (Figure 4A). The highest chromatin accessibility for RCF-specific loci was observed in CD200⁺/CD146⁻ CFs. Furthermore, a motif enrichment analysis identified several TFs such as RUNX1, SMAD, AP-1(JUN), and TEAD as candidate regulators of RCFs (Figure 4B). With the exception of Runx1, most of the single-cell expression patterns of the identified TFs showed no specificity for RCFs (Figure XIIA in the Data Supplement). Although Runx1 has been reported to play a role in CF activation,³¹ overexpression of Runx1 in cultured CFs was unable to induce the RCFs-specific transcription signature (Figure 4C). As a result, we considered that the analysis and the characterization of the TF role should be extended to compare the ATAC-seq fingerprint between CD200⁺/CD146⁻ and CD200⁻/CD146⁻ (Figure XIIB in the Data Supplement). Based on motifs detected in peaks from both populations and in the percentage of increased or decreased accessibility, TFs were grouped in 3 different sets: enriched in CD200⁺ peaks but not differentially expressed (set 1); enriched in CD200⁻ peaks but not differentially expressed (set 2); and enriched in CD200⁻ peaks and differentially expressed (set 3). No TFs enriched in CD200⁺ peaks and differentially expressed were found (Figure XIIC in the Data Supplement). To characterize TFs in sets 1 and 3, we used peak-gene mapping to conduct a gene-set analysis per TF. These analyses revealed specific TFs in set 1 (ATF3, JUN, or ZNF93) and in set 3 (RUNX1, WT1, or KLF5; Figure XIID in the Data Supplement). Together, these results represent an initial approach for the characterization of the TFs involved in the generation of RCFs.

• Download figure

• Download PowerPoint

Figure 4. Molecular regulation of reparative cardiac fibroblasts (RCFs) identity. A, Genome browser snapshots show accessibility profiles of representative loci in the global cardiac fibroblast (CF) population at different time points and in CD200⁺ and CD200⁻ CF subpopulations at 7 days postinfarction (dpi). Shadowed areas mark distal regulatory elements displaying increased accessibility in the CD200⁺ subpopulation. B, Dot plot represents motif enrichment and expression specificity of potential transcription factor (TF)–mediating RCF response (cluster B/CD200⁺ CF). First row, analysis of CD200⁺-specific accessible distal regulatory elements (±1.5 Kb from the transcriptional start site); second row, analysis of distal regulatory elements found within RCF-specific loci; third row, expression values of TF in single-cell RNA sequencing (scRNA-seq); fourth row, expression values for TF in CD200⁺ bulk RNA sequencing. C, Volcano plot shows differential gene expression of in vitro grown CF overexpressing Runx1. Genes with log fold change of ±1.5 and P<0.05 were considered differentially expressed. Red dots represent RCF markers. D, Transforming growth factor (TGF)–β network centrality analysis revealed noncanonical PI3K-Akt pathway related with RCF markers. E, Heatmap shows relative expression of RCF markers in nontreated (control), TGF-β, TGF-β + vehicle (dimethyl sulfoxide [DMSO]), and TGF-β + LY294002 cultured CFs. Linear model differential expression, n=4 per group. F, Quantification of area covered by cultured CFs after 23 hours of treatment in wound-healing experiment. One-way analysis of variance, Kruskal-Wallis post hoc test, n=4 to 8. **P≤0.01, ***P≤0.001. ATAC-seq indicates assay for transposase-accessible chromatin sequencing; and GFP, green fluorescent protein.

As an additional effort to characterize molecular regulation of RCFs, we built a network spanning protein signaling and gene-regulatory interactions using our single-cell expression dataset based on the distances between cluster markers and signaling receptors. This analysis highlighted the noncanonical TGF-β1/PI3K-Akt pathway over the canonical one as the main driver of RCF gene expression (Figure 4D). To validate this, we stimulated in vitro CFs with TGF-β in the presence or absence of LY294002, a known PI3K-Akt inhibitor.³² This inhibitor abolished the transcriptomic signature of RCFs induced by TGF-β (Figure 4E and Figure XIIIA and XIIIB in the Data Supplement) as well as the ability of fibroblasts to migrate and proliferate (Figure 4F and the Movie I in the Data Supplement). These results emphasize the relevance of the noncanonical TGF-β1/PI3K-Akt pathway in controlling identity and function of RCFs.

To explore the role of TGF-β in depth, we cultured adult CFs in the presence of small molecules that specifically inhibit different elements of the TGF-β1 signaling pathway (Table II in the Data Supplement and Figure XIIIC in the Data Supplement). The inhibition of specific elements of the canonical (SMAD2/3)^33,34 and noncanonical (AKT and p38)^33,35–38 TGF-β1 signaling pathways reduced the expression of both ECM (Col1α1, Col3α1, Lox) and CF activation (Cthrc1, Ddah1, Postn, Acta2)–related genes (Figure XIIID in the Data Supplement). However, PD98059, a specific inhibitor of ERK,^33,34,38 did not affect any of the genes (Figure XIIID in the Data Supplement). These results indicate that activation of the profibrotic profile and the phenotypic transformation of RCFs are regulated by a balance between SMAD2/3 and AKT/p38, but not via ERK. Remarkably, SB-431542, a specific inhibitor of TGF-β receptor 1/ALK5,^33,38,39 selectively reduced the RCF markers together with Acta2 and Col1α1, but not Postn and Col3α1. This highlights the complexity that underlies the regulation of CF activation.

CTHRC1-Mediated RCF Activity Is Essential for the Formation of the Healing Scar

As Cthrc1 is exclusively expressed in CFs but not in other cardiac cell populations such as CD45, CD31, or cardiomyocytes (Figure 5A through 5C and Figure XIVB and XIVC in the Data Supplement), we decided to analyze the functional role of RCFs in mice with genetic ablation of Cthrc1 (Cthrc1-KO; Figure XIVA in the Data Supplement). Cthrc1-KO mice appear normal during development and adulthood. No changes were observed in size or weight, and more importantly, no cardiac phenotype was detected.²² However, after MI, knock-out (KO) mice showed a dramatic decrease in survival (80% wild-type [WT] versus 30% KO) attributable to ventricular rupture (Figure 5D). This was associated with a significant decrease in collagen deposition in the free wall of the left ventricle of KO infarcted hearts (50%; Figure 5E). This phenotype was not associated with differences in the transcriptomic profile of other cardiac cell populations between WT and KO mice, suggesting that CTHRC1 had no pleiotropic effects (Figure XIVC in the Data Supplement). In contrast, we found significant differences between KO and WT CFs in gene ontology terms related to cell division, proliferation, and protein synthesis on MI (Figure 5F and 5G). Cultured KO CFs from healthy mice stimulated with TGF-β1 showed a decrease in gene ontology terms related to angiogenesis, muscle contraction, and vasculature development when compared with WT CFs, indicating a reduced capacity to respond to TGF-β1 (Figure XV in the Data Supplement).

• Download figure

• Download PowerPoint

Figure 5. CTHRC1 is an essential effector of reparative cardiac fibroblasts (RCFs) for the healing repair process. A, Normalized expression bar plot (mean±SD) of Cthrc1 in cardiac fibroblasts (CFs; GFP⁺/CD31⁻/CD45⁻), endothelial cells (CD31⁺), and bone marrow–derived cells (CD45⁺) at different time points. B, Localization of CTHRC1 (green) in the left ventricle of healthy wild-type (WT) mice and at 3, 5, and 7 days postinfarction (dpi). Cardiac troponin-I⁺ cardiomyocytes, red; nuclei (DAPI [4′,6-diamidino-2-phenylindole]), blue. C, Normalized expression bar plot (mean±SD) of Cthrc1 in endothelial cells (CD31⁺), bone marrow–derived cells (CD45⁺), CFs (mEFSK4⁺/CD31⁻/CD45⁻), and cardiomyocytes from WT hearts at 5 dpi; linear model with normalized counts. D, Kaplan-Meier survival curves after myocardial infarction (MI) in WT and Cthrc1 knockout (KO) mice. Log-rank Mantel-Cox test, n=10 per group. E, Representative images of collagen deposition in left ventricle of WT (left) and KO (right) mice. Quantification of collagen deposition in the left ventricle (LV) in both genotypes at 3 dpi (WT, open circles; KO, closed circles), unpaired t test (right). F, Volcano plot shows differential gene expression between WT and KO CFs at 5 dpi. G, Dot plot comparison of enriched pathways in the bulk RNA sequencing analysis between Cthrc1-KO (left column) and WT (right column) at 7 dpi. H, t-Distributed stochastic neighbor embedding (t-SNE) representation of 4189 CFs from 1 KO heart at 7 dpi. Red dots represent RCF-like fibroblasts and the number the percentage of them in relation to the total isolated CF. t-SNE representation for RCF markers (below). I, Proportion of cluster B CFs in each of the datasets. **P≤0.01, ***P≤0.001. SRP indicates signal recognition particle.

scRNA-seq analysis on Cthrc1-KO hearts at 7 dpi (4189 cells) did not show a reduction of RCFs but rather an increased percentage of RCF-like cells in comparison with WT mice (21% versus 14%; Figure 5H and 5I). RCF-like CFs and other subpopulations of activated CFs such as POSTN⁺ were found in the IZ and BZ of Cthrc1-KO infarcted mice, indicating that the ventricular rupture phenotype is derived from the absence of CTHRC1 (Figure XIVD and XIVE in the Data Supplement). Accordingly, deletion of Cthrc1 was associated with downregulation of genes related to ECM organization, collagen biosynthesis, and TGF-β–mediated regulation of the ECM in CFs in comparison with WT mice (Figure XIVF in the Data Supplement). Globally, these results suggested that RCFs orchestrate cardiac repair via secretion of CTHRC1, which affects the deposition and synthesis of ECM molecules and promotes the proliferation of CFs.

An RCF-Like Expression Signature Is Detected in a Preclinical Model of MI and in Human Cardiac Fibrotic Tissue

To assess the translational potential of our findings, we examined whether the RCFs signature can be detected in a preclinical pig model of MI (Figure 6A). This was done by comparing the global transcriptome of representative biopsies from IZ and remote zone of infarcted pig hearts at 8, 60, and 180 dpi (Figure 6B). Despite the lack of cell type resolution, we were able to detect the induction of 33%, 42%, and 32% of the RCF signature at 8, 60, and 180 dpi, respectively, in the IZ but not in the remote zone (Figure XVIA in the Data Supplement). Consistent with murine data, CTHRC1⁺ cells localized to the IZ but not to the remote zone, in the same region where deposits of COL1α1 and POSTN were found (Figure 6C). CTHRC1 expression was lower at 180 dpi in comparison with 8 and 60 dpi (Figure 6B) with no CTHRC1⁺ cells detected at 180 dpi by immunohistochemistry (data not shown). At 8 dpi, we identified 3 different subpopulations of activated CFs—CTHRC1⁺, αSMA⁺, and αSMA⁺/CTHRC1⁺—in the IZ and BZ, with a decrease at 60 dpi, corroborating the heterogeneity of activated CFs observed in mice (Figure XVIB in the Data Supplement). These results indicated that an RCFs-like population also appears after MI in our preclinical swine model with similar dynamics as in mice.

• Download figure

• Download PowerPoint

Figure 6. Reparative cardiac fibroblasts (RCFs) marker expression correlates with cardiac function in infarcted pigs and is conserved in humans. A, Twenty-nine pigs underwent ischemia/reperfusion surgery and 2 animals were used as controls. Cardiac function was determined in 21 pigs after 6 months using echocardiography and MRI. B, Normalized expression bar plots (mean±SD) for zonal transcriptomic profiling at 8 days postinfarction (dpi; top; n=2), 60 dpi (middle; n=6), and 180 dpi (bottom; n=9). Normalized expression bar plots (mean±SD) for top RCF markers for healthy hearts (blue) and hearts at 8, 60, and 180 dpi in remote zone (RZ; green) and infarct zone (IZ; red). C, Immunohistochemistry of COL1α1, POSTN, and CTHRC1 in IZ and RZ at 8 dpi. D, Distribution of pigs with reduced (≤45%, circles; n=9) or preserved (>45%, squares; n=12) ejection fraction (EF) at 180 dpi by echocardiography or MRI. Two-tailed t test with a Mann-Whitney post hoc test. E, Comparison of POSTN, COL1a1, and CTHRC1 expression levels in different anatomic regions between infarcted pigs with reduced or preserved EF. Two-way analysis of variance, n=8 to 10 (top). Correlation between EF and expression level of POSTN, COL1a1, and CTHRC1 in IZ. Spearman correlation coefficient (bottom). F, Distribution of pigs with reduced (≤45%, circles) or preserved (>45%, squares) infarct area at 180 dpi by magnetic resonance imaging. Two-tailed t test with a Mann-Whitney post hoc test, n=21. G, Correlations between infarct area and expression level of POSTN, COL1a1, and CTHRC1 in IZ. Spearman correlation coefficient. H, Normalized expression bar plots for top RCF markers (CTHRC1, DDAH1, POSTN, FMOD, LOX, PTN, and COMP) in human samples from healthy hearts (left ventricle [LV] and right ventricle [RV]; n=6), infarcted hearts (IZ and RZ; n=8), and hearts with dilated cardiomyopathy (DCM; LV and RV; n=5). Likelihood ratio test. Principal component (PC) analysis scatterplot of human samples subjected to transcriptomic profiling. I, Representative images of the immunohistochemistry analysis of CTHRC1⁺ CFs (in brown; arrows) performed on sections of the IZ and RZ obtained from the LV of 2 different patients who had myocardial infarction (MI). BZ indicates border zone. *P≤0.05, **P≤0.01, ***P≤0.001, ****P≤0.0001. Scale bars, 100 µm.

To assess the relation between the expression of CTHRC1 and cardiac function, we compared the level of expression of CTHRC1 in the IZ between animals with midrange, preserved (>45%), and reduced (≤45%) ejection fraction⁴⁰ at 180 dpi (n=12 versus n=9, respectively; Figure 6D). Despite the critical role of CTHRC1 in the early stages after MI, we found that at later stages (180 dpi) the level of CTHRC1 (BZ, IZ) was significantly higher in the group with cardiac dysfunction (Figure 6E). Consequently, a negative correlation between the ejection fraction and the expression of CTHRC1, COL1α1, and POSTN in the IZ was found (Figure 6E), which was only statistically significant for CTHRC1 and COL1α1. A positive and significant correlation between the infarct area and expression levels was observed for CTHRC1 and POSTN but not for COL1α1 (Figure 6F and 6G). Taken together, our results suggest that the upregulation of CTHRC1 could be beneficial immediately after MI, whereas increased expression at later time points may reflect incomplete repair with progressive fibrosis.

The transcriptomic signature of RCFs was partially found in biopsies obtained from the ischemic zone of patients with ischemic cardiomyopathy and in biopsies from both right and left ventricle of patients with dilated cardiomyopathy. Several genes of the RCF transcriptomic signature (CTHRC1, PTM, FMOD) were significantly overexpressed in all the pathologic conditions included in this study compared with controls, as shown in the bar plots from Figure 6H. Moreover, CTHRC1⁺ CFs were found in the ischemic zone but not in remote zone of infarcted hearts from patients (Figure 6I). These findings highlight the potential role of CTHRC1⁺ RCFs in orchestrating the cardiac repair process in patients with cardiac fibrosis.

Discussion

As the main cellular component responsible for cardiac healing after MI, CFs are an attractive therapeutic target. However, the lack of appropriate markers, the cellular heterogeneity, and the limited understanding of the molecular mechanisms underlying their activation has precluded the development of successful therapies targeting fibrotic remodeling.⁴ In this study, we characterize the heterogeneity of CFs using Col1α1-GFP mice at the single-cell level. We identify a specific subpopulation of CFs defined by the expression of a unique transcriptional signature responsible for the healing response after cardiac injury. Among the top markers, we identified Cthrc1, a crucial molecule involved in the ventricular remodeling process. In a recent study, a similar role for Cthrc1⁺ fibroblasts has been described in fibrotic lungs.⁴¹ We also identified RCFs in a large swine model and in patients with cardiac diseases that lead to myocardial fibrosis.

Nonmyocyte cardiac cell heterogeneity has been recently characterized with a particular focus on CFs.^7–10,12 Most of these studies focus on the early phase after cardiac injury¹² or assess other pathologic conditions.^7,9,10 In the study by Kretzschmar and colleagues,¹⁰ the transcriptome of 282 CFs at 14 dpi was studied, identifying 11 clusters, 1 of which was defined as activated CFs (FstlI⁺). This subpopulation had similar properties as RCFs reported in this study but both subpopulations showed different transcriptomic signatures. Fu et al⁸ recently described a new subpopulation of activated CFs, the matrifibrocytes, characterized by the expression of ECM and tendon genes. Despite expression of several RCF marker genes in the matrifibrocytes, it was not possible to compare both subpopulations at the single-cell level. Moreover, RCFs have a role during the healing scar, in contrast to matrifibrocytes, with their role in supporting the mature scar at later stages. Although markers associated with RCFs have been described in other CF subpopulations,^7–9,12 our transcriptomic, histologic, and functional results clearly indicate that RCFs are distinct from Postn⁺, Acta2⁺, or matrifibrocytes. An important difference between our work and previous studies is that we focused exclusively on the CF population, unlike recently reported scRNA-seq studies in models of MI.^9,10,12 This distinction together with the careful description of the anatomic location and the temporal dynamics of RCFs after MI make for a better characterization of CFs heterogeneity and a refined definition of new subpopulations.

Understanding the regulation of RCFs will be essential for the identification of therapeutic strategies. On the one hand, our transcriptomic studies together with our binding site analyses indicate a potential role for different TFs in the generation of RCFs. We validated the role of SOX9, consistent with recent publications.^28,42 Moreover, our ATAC-seq data revealed TFs that could be involved in the regulation of RCF gene markers and identify other TFs involved in the regulation of CD200⁺ cells through chromatin accessibility.

Our in silico approaches emphasize the role of the TGF-β signaling pathway as a master regulator of CF activation and the development of cardiac fibrosis.^43,44 TGF-β signaling is known to promote myofibroblast formation and ECM production.^{28,29,45–48} Our results demonstrate that both canonical and noncanonical TGF-β signaling pathways have a significant contribution in the activation of RCFs, with a predominant role for the PI3K-Akt pathway in generating RCFs. This pathway upregulates Cthrc1, the key marker gene for RCF, and the migratory and proliferative capacities of CFs. These results suggest that the development of specific inhibitors of the PI3K-Akt pathway could define new targets for modulating cardiac fibrosis, although further studies of the downstream mechanisms are required.

CTHRC1 is a secreted protein that participates in collagen matrix synthesis through TGF-β signaling.^21,22,49 Moreover, the expression of Cthrc1 has been recently identified as a potential marker for activated CFs in the heart^7,9,12 and lungs.⁴¹ The longer follow-up in our study and the evaluation of Cthrc1 deletion in the context of MI may explain why it was not previously identified as a key player in the cardiac repair process. Consistent with other studies where key components of activated CFs were ablated,^7,10 the deletion of Cthrc1 was associated with increased mortality attributable to ventricular rupture. In our case, RCFs were still detected in the myocardium of KO mice, together with other subpopulations of activated CFs. Moreover, our in silico data indicated that RCFs originate from activated Postn⁺ CFs with a specific role in the secretion and deposition of collagen in the healing scar. These results suggest that the ventricular rupture described by others could be an effect of the deletion of CFs that give rise to RCFs.^36,50,51 Taken together, these data indicate that Cthrc1 may be an effector of RCFs, promoting the healing process, rather than an inducer of the activation of this population.

Our results in a preclinical model of MI and in myocardial biopsy specimens of patients with ischemic or dilated cardiomyopathy provide additional translational relevance to our findings. Similar to the observations in mice, expression of CTHRC1 was upregulated in the fibrotic areas with identification of CTHRC1⁺ CFs in these zones and not in the healthy myocardium. Also, the pattern of expression after MI in pigs showed a significant increase in the level of CTHRC1 in the early stages after MI, followed by a reduction of its expression levels in the chronic phase. Of note, the presence of more CTHRC1⁺ CFs or higher levels of CTHRC1 expression, particularly at later time points, was associated with worse cardiac function, which would seem counterintuitive because CTHRC1 expression is associated with cardiac healing. These findings suggest that both responses depend on ECM deposition by the CTHRC1⁺ CF lineage, such that ECM contributes to acute cardiac repair and prevents rupture, whereas sustained ECM deposition at later stages may reflect pathologic remodeling. Although additional studies are needed, our findings suggest this molecule might be considered as a potential candidate to study “adverse remodeling” versus “myocardial recovery” in heart failure, based on 3 properties of the molecule: the dynamic transcriptome of CTHRC1; its essential role on ECM during the early healing response; and its potentially negative effect on cardiac function during later stages of adverse ventricular remodeling.⁵²

Our study, including the largest number of CFs interrogated by scRNA-seq, clearly demonstrates the remarkable heterogeneity of CFs after MI and defines a specific subpopulation, RCFs, characterized by the expression of Cthrc1. Insights into molecular regulation and biological function distinctly support the role of Cthrc1⁺ CFs in the early healing process after MI. These novel findings should facilitate the search for more personalized treatment through control of the cardiac fibrotic process.

Study Limitations

The current analysis does not provide complete mechanistic insight because regulation based on posttranslational modifications of TF was not considered. With the apparent complexity underlying the regulation of CFs, we cannot rule out that other molecules or subpopulations of CF could have an essential role in the ventricular remodeling process. Future studies should address whether deletion of RCFs using Cre-inducible diphtheria toxin receptor transgenic mouse models instead of CTHRC1 null mice results in a different outcome.⁵³

Competency in Medical Knowledge

Unraveling CFs heterogeneity for understanding different cardiac fibrotic processes is critical. In the case of MI, Cthrc1⁺ CFs are responsible for a proper healing response to cardiac injury. Its role in chronic models requires further studies.

Translational Outlook

Further investigations are needed to determine the potential role of CTHRC1 during MI related to cardiac repair and its mechanism of regulation to identify the therapeutic windows that may allow for controlling the size of the fibrotic scar.

Acknowledgments

The authors thank D.A. Brenner and T. Kisseleva (University of California San Diego) for the gift of Col1α1-GFP mice; S. Sarvide for her comments in the preparation of this article; and MMC BioBank, a core facility of Maine Medical Center Research Institute.

Sources of Funding

This work was supported by Instituto de Salud Carlos III and Fondo Europeo de Desarrollo Regional funds (PI16/00129, CPII15/00017, PI19/00501), Red de Terapia Celular RD16/0011/0005 and Ministerio de Economía y Empresa (Program RETOS Cardiomesh), ERANET II (Nanoreheart), and the Horizon 2020 Program BRAVE. Dr Ruiz-Villalba is supported by Fondo Social Europeo/Ministerio de Economía, Industria y Competitividad–Agencia Estatal de Investigación/IJCI-2016-30254, and the Spanish Ministerio de Ciencia, Innovación y Universidades (RTI2018-095410-BI00). Dr Fortelny is supported by a fellowship from the European Molecular Biology Organization (EMBO ALTF 241–2017). Dr Bock is supported by a New Frontiers Group award of the Austrian Academy of Sciences and by a European Research Council Starting Grant (EU Horizon 2020 research and innovation program, grant agreement 679146). This work was supported in part by the National Heart, Lung, and Blood Institute of the National Institutes of Health under grants R01 HL136560 and HL146504 (Dr Lindner and Dr Ryzhov). This work also used services of the Physiology Core, Histopathology Core, Confocal Microscopy Core, and Mouse Transgenic Core, which are supported by National Institutes of Health/National Institute of General Medical Sciences grants P30GM106391 and P20GM121301.

Disclosures

None.

Supplemental Materials

Methods

Data Supplement Tables I–VI

Data Supplement Figures I–XVI

Data Supplement Movie I

File

Data Resource: GSE132146

References 54–89

Footnotes

*Drs Ruiz-Villalba, Romero, and Hernandez contributed equally and should be considered joint first authors.

†Drs Lindner, Lara-Astiaso, and Prósper contributed equally and should be considered joint senior authors.

Sources of Funding, see page 1845

https://www./journal/circ

The Data Supplement is available with this article at https://www./doi/suppl/10.1161/circulationaha.119.044557.

Felipe Prosper, MD, Department of Hematology and Cell Therapy, Clínica Universidad de Navarra, Pamplona, 31008, Spain. Email [e-mail]

David Lara-Astiaso, PhD, Advanced Genomics Laboratory, Program of Hemato-Oncology, CIMA, University of Navarra, Pamplona, 31008, Spain. Email [e-mail]

Volkhard Lindner, MD, PhD, Maine Medical Centre Research Institute, Scarborough, ME 04072. Email [e-mail]

References

• 1.

  Pinto AR, Ilinykh A, Ivey MJ, Kuwabara JT, D’Antoni ML, Debuque R, Chandran A, Wang L, Arora K, Rosenthal NA, et al.

  . Revisiting cardiac cellular composition.Circ Res. 2016; 118:400–409. doi: 10.1161/CIRCRESAHA.115.307778LinkGoogle Scholar

• 2.

  Frangogiannis NG

  . Pathophysiology of myocardial infarction.Compr Physiol. 2015; 5:1841–1875. doi: 10.1002/cphy.c150006CrossrefMedlineGoogle Scholar

• 3.

  Shinde AV, Frangogiannis NG

  . Mechanisms of fibroblast activation in the remodeling myocardium.Curr Pathobiol Rep. 2017; 5:145–152. doi: 10.1007/s40139-017-0132-zCrossrefMedlineGoogle Scholar

• 4.

  Tallquist MD, Molkentin JD

  . Redefining the identity of cardiac fibroblasts.Nat Rev Cardiol. 2017; 14:484–491. doi: 10.1038/nrcardio.2017.57CrossrefMedlineGoogle Scholar

• 5.

  Ruiz-Villalba A, Simón AM, Pogontke C, Castillo MI, Abizanda G, Pelacho B, Sánchez-Domínguez R, Segovia JC, Prósper F, Pérez-Pomares JM

  . Interacting resident epicardium-derived fibroblasts and recruited bone marrow cells form myocardial infarction scar.J Am Coll Cardiol. 2015; 65:2057–2066. doi: 10.1016/j.jacc.2015.03.520CrossrefMedlineGoogle Scholar

• 6.

  Moore-Morris T, Guimarães-Camboa N, Banerjee I, Zambon AC, Kisseleva T, Velayoudon A, Stallcup WB, Gu Y, Dalton ND, Cedenilla M, et al.

  . Resident fibroblast lineages mediate pressure overload-induced cardiac fibrosis.J Clin Invest. 2014; 124:2921–2934. doi: 10.1172/JCI74783CrossrefMedlineGoogle Scholar

• 7.

  Kanisicak O, Khalil H, Ivey MJ, Karch J, Maliken BD, Correll RN, Brody MJ, Lin S-CJ, Aronow BJ, Tallquist MD, et al.

  . Genetic lineage tracing defines myofibroblast origin and function in the injured heart.Nat Commun. 2016; 7:12260. doi: 10.1038/ncomms12260CrossrefMedlineGoogle Scholar

• 8.

  Fu X, Khalil H, Kanisicak O, Boyer JG, Vagnozzi RJ, Maliken BD, Sargent MA, Prasad V, Valiente-Alandi I, Blaxall BC, et al.

  . Specialized fibroblast differentiated states underlie scar formation in the infarcted mouse heart.J Clin Invest. 2018; 128:2127–2143. doi: 10.1172/JCI98215CrossrefMedlineGoogle Scholar

• 9.

  Gladka MM, Molenaar B, de Ruiter H, van der Elst S, Tsui H, Versteeg D, Lacraz GPA, Huibers MMH, van Oudenaarden A, van Rooij E

  . Single-cell sequencing of the healthy and diseased heart reveals cytoskeleton-associated protein 4 as a new modulator of fibroblasts activation.Circulation. 2018; 138:166–180. doi: 10.1161/CIRCULATIONAHA.117.030742LinkGoogle Scholar

• 10.

  Kretzschmar K, Post Y, Bannier-Hélaouët M, Mattiotti A, Drost J, Basak O, Li VSW, van den Born M, Gunst QD, Versteeg D, et al.

  . Profiling proliferative cells and their progeny in damaged murine hearts.Proc Natl Acad Sci U S A. 2018; 115:E12245–E12254. doi: 10.1073/pnas.1805829115CrossrefMedlineGoogle Scholar

• 11.

  Skelly DA, Squiers GT, McLellan MA, Bolisetty MT, Robson P, Rosenthal NA, Pinto AR

  . Single-cell transcriptional profiling reveals cellular diversity and intercommunication in the mouse heart.Cell Rep. 2018; 22:600–610. doi: 10.1016/j.celrep.2017.12.072CrossrefMedlineGoogle Scholar

• 12.

  Farbehi N, Patrick R, Dorison A, Xaymardan M, Janbandhu V, Wystub-Lis K, Ho JW, Nordon RE, Harvey RP

  . Single-cell expression profiling reveals dynamic flux of cardiac stromal, vascular and immune cells in health and injury.Elife2019; 8:e43882. doi: 10.7554/eLife.43882CrossrefMedlineGoogle Scholar

• 13.

  Pyagay P, Heroult M, Wang Q, Lehnert W, Belden J, Liaw L, Friesel RE, Lindner V

  . Collagen triple helix repeat containing 1, a novel secreted protein in injured and diseased arteries, inhibits collagen expression and promotes cell migration.Circ Res. 2005; 96:261–268. doi: 10.1161/01.RES.0000154262.07264.12LinkGoogle Scholar

• 14.

  Jin YR, Stohn JP, Wang Q, Nagano K, Baron R, Bouxsein ML, Rosen CJ, Adarichev VA, Lindner V

  . Inhibition of osteoclast differentiation and collagen antibody-induced arthritis by CTHRC1.Bone. 2017; 97:153–167. doi: 10.1016/j.bone.2017.01.022CrossrefMedlineGoogle Scholar

• 15.

  Yata Y, Scanga A, Gillan A, Yang L, Reif S, Breindl M, Brenner DA, Rippe RA

  . DNase I-hypersensitive sites enhance alpha1(I) collagen gene expression in hepatic stellate cells.Hepatology. 2003; 37:267–276. doi: 10.1053/jhep.2003.50067CrossrefMedlineGoogle Scholar

• 16.

  Ivey MJ, Kuwabara JT, Pai JT, Moore RE, Sun Z, Tallquist MD

  . Resident fibroblast expansion during cardiac growth and remodeling.J Mol Cell Cardiol. 2018; 114:161–174. doi: 10.1016/j.yjmcc.2017.11.012CrossrefMedlineGoogle Scholar

• 17.

  Ivey MJ, Tallquist MD

  . Defining the cardiac fibroblast.Circ J. 2016; 80:2269–2276. doi: 10.1253/circj.CJ-16-1003CrossrefMedlineGoogle Scholar

• 18.

  Horn MA, Trafford AW

  . Aging and the cardiac collagen matrix: Novel mediators of fibrotic remodelling.J Mol Cell Cardiol. 2016; 93:175–185. doi: 10.1016/j.yjmcc.2015.11.005CrossrefMedlineGoogle Scholar

• 19.

  Gonzalez A, Schelbert EB, Diez J, Butler J

  . Myocardial interstitial fibrosis in heart failure: biological and translational perspectives.J Am Coll Cardiol2018; 71:1696–1706. doi: 10.1016/j.jacc.2018.02.021CrossrefMedlineGoogle Scholar

• 20.

  Binks AP, Beyer M, Miller R, LeClair RJ

  . Cthrc1 lowers pulmonary collagen associated with bleomycin-induced fibrosis and protects lung function.Physiol Rep. 2017; 5:e13115. doi: 10.14814/phy2.13115CrossrefMedlineGoogle Scholar

• 21.

  LeClair RJ, Durmus T, Wang Q, Pyagay P, Terzic A, Lindner V

  . Cthrc1 is a novel inhibitor of transforming growth factor-beta signaling and neointimal lesion formation.Circ Res. 2007; 100:826–833. doi: 10.1161/01.RES.0000260806.99307.72LinkGoogle Scholar

• 22.

  Stohn JP, Perreault NG, Wang Q, Liaw L, Lindner V

  . Cthrc1, a novel circulating hormone regulating metabolism.PLoS One. 2012; 7:e47142. doi: 10.1371/journal.pone.0047142CrossrefMedlineGoogle Scholar

• 23.

  Wang Y, Lee M, Yu G, Lee H, Han X, Kim D

  . CTHRC1 activates pro-tumorigenic signaling pathways in hepatocellular carcinoma.Oncotarget2017; 8:105238. doi: 10.18632/oncotarget.22164CrossrefMedlineGoogle Scholar

• 24.

  Xu G, Fan W, Wang F, Lu H, Xing X, Zhang R, Jiang P

  . CTHRC1 as a novel biomarker in the diagnosis of cervical squamous cell carcinoma.Int J Clin Exp Pathol. 2018; 11:847–854.MedlineGoogle Scholar

• 25.

  La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, Lidschreiber K, Kastriti ME, Lönnerberg P, Furlan A, et al.

  . RNA velocity of single cells.Nature. 2018; 560:494–498. doi: 10.1038/s41586-018-0414-6CrossrefMedlineGoogle Scholar

• 26.

  Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ

  . Generalizing RNA velocity to transient cell states through dynamical modeling [published online August 3, 2020].Nat Biotechnol. doi: 10.1038/s41587-020-0591-3. https://www./articles/s41587-020-0591-3Google Scholar

• 27.

  Hu HH, Chen DQ, Wang YN, Feng YL, Cao G, Vaziri ND, Zhao YY

  . New insights into TGF-beta/Smad signaling in tissue fibrosis.Chem Biol Interact2018; 292:76–83. doi: 10.1016/j.cbi.2018.07.008CrossrefMedlineGoogle Scholar

• 28.

  Lacraz GPA, Junker JP, Gladka MM, Molenaar B, Scholman KT, Vigil-Garcia M, Versteeg D, de Ruiter H, Vermunt MW, Creyghton MP, et al.

  . Tomo-Seq identifies SOX9 as a key regulator of cardiac fibrosis during ischemic injury.Circulation. 2017; 136:1396–1409. doi: 10.1161/CIRCULATIONAHA.117.027832LinkGoogle Scholar

• 29.

  Khalil H, Kanisicak O, Prasad V, Correll RN, Fu X, Schips T, Vagnozzi RJ, Liu R, Huynh T, Lee SJ, et al.

  . Fibroblast-specific TGF-beta-Smad2/3 signaling underlies cardiac fibrosis.J Clin Invest2017; 127:3770–3783. doi: 10.1172/JCI94753CrossrefMedlineGoogle Scholar

• 30.

  Huebener P, Abou-Khamis T, Zymek P, Bujak M, Ying X, Chatila K, Haudek S, Thakker G, Frangogiannis NG

  . CD44 is critically involved in infarct healing by regulating the inflammatory and fibrotic response.J Immunol. 2008; 180:2625–2633. doi: 10.4049/jimmunol.180.4.2625CrossrefMedlineGoogle Scholar

• 31.

  Kim W, Barron DA, San Martin R, Chan KS, Tran LL, Yang F, Ressler SJ, Rowley DR

  . RUNX1 is essential for mesenchymal stem cell proliferation and myofibroblast differentiation.Proc Natl Acad Sci U S A. 2014; 111:16389–16394. doi: 10.1073/pnas.1407097111CrossrefMedlineGoogle Scholar

• 32.

  Montiel-Duarte C, Cordeu L, Agirre X, Román-Gómez J, Jiménez-Velasco A, José-Eneriz ES, Gárate L, Andreu EJ, Calasanz MJ, Heiniger A, et al.

  . Resistance to imatinib mesylate-induced apoptosis in acute lymphoblastic leukemia is associated with PTEN down-regulation due to promoter hypermethylation.Leuk Res. 2008; 32:709–716. doi: 10.1016/j.leukres.2007.09.005CrossrefMedlineGoogle Scholar

• 33.

  Vivar R, Humeres C, Ayala P, Olmedo I, Catalán M, García L, Lavandero S, Díaz-Araya G

  . TGF-β1 prevents simulated ischemia/reperfusion-induced cardiac fibroblast apoptosis by activation of both canonical and non-canonical signaling pathways.Biochim Biophys Acta. 2013; 1832:754–762. doi: 10.1016/j.bbadis.2013.02.004CrossrefMedlineGoogle Scholar

• 34.

  Voloshenyuk TG, Landesman ES, Khoutorova E, Hart AD, Gardner JD

  . Induction of cardiac fibroblast lysyl oxidase by TGF-β1 requires PI3K/Akt, Smad3, and MAPK signaling.Cytokine. 2011; 55:90–97. doi: 10.1016/j.cyto.2011.03.024CrossrefMedlineGoogle Scholar

• 35.

  Ellis IR, Jones SJ, Lindsay Y, Ohe G, Schor AM, Schor SL, Leslie NR

  . Migration stimulating factor (MSF) promotes fibroblast migration by inhibiting AKT.Cell Signal. 2010; 22:1655–1659. doi: 10.1016/j.cellsig.2010.06.005CrossrefMedlineGoogle Scholar

• 36.

  Molkentin JD, Bugg D, Ghearing N, Dorn LE, Kim P, Sargent MA, Gunaje J, Otsu K, Davis J

  . Fibroblast-specific genetic manipulation of p38 mitogen-activated protein kinase in vivo reveals its central regulatory role in fibrosis.Circulation. 2017; 136:549–561. doi: 10.1161/CIRCULATIONAHA.116.026238LinkGoogle Scholar

• 37.

  Papakrivopoulou J, Lindahl GE, Bishop JE, Laurent GJ

  . Differential roles of extracellular signal-regulated kinase 1/2 and p38MAPK in mechanical load-induced procollagen α1 (I) gene expression in cardiac fibroblasts.Cardiovasc Res2004; 61:736–744. doi: 10.1016/j.cardiores.2003.12.018CrossrefMedlineGoogle Scholar

• 38.

  Hu J, Wang X, Wei SM, Tang YH, Zhou Q, Huang CX

  . Activin A stimulates the proliferation and differentiation of cardiac fibroblasts via the ERK1/2 and p38-MAPK pathways.Eur J Pharmacol. 2016; 789:319–327. doi: 10.1016/j.ejphar.2016.07.053CrossrefMedlineGoogle Scholar

• 39.

  Bollong MJ, Yang B, Vergani N, Beyer BA, Chin EN, Zambaldo C, Wang D, Chatterjee AK, Lairson LL, Schultz PG

  . Small molecule-mediated inhibition of myofibroblast transdifferentiation for the treatment of fibrosis.Proc Natl Acad Sci U S A. 2017; 114:4679–4684. doi: 10.1073/pnas.1702750114CrossrefMedlineGoogle Scholar

• 40.

  Ponikowski P, Voors AA, Anker SD, Bueno H, Cleland JGF, Coats AJS, Falk V, González-Juanatey JR, Harjola VP, Jankowska EA, et al.

  ; ESC Scientific Document Group. 2016 ESC Guidelines for the diagnosis and treatment of acute and chronic heart failure: the task force for the diagnosis and treatment of acute and chronic heart failure of the European Society of Cardiology (ESC): developed with the special contribution of the Heart Failure Association (HFA) of the ESC.Eur Heart J. 2016; 37:2129–2200. doi: 10.1093/eurheartj/ehw128CrossrefMedlineGoogle Scholar

• 41.

  Tsukui T, Sun KH, Wetter JB, Wilson-Kanamori JR, Hazelwood LA, Henderson NC, Adams TS, Schupp JC, Poli SD, Rosas IO, et al.

  . Collagen-producing lung cell atlas identifies multiple subsets with distinct localization and relevance to fibrosis.Nat Commun. 2020; 11:1920. doi: 10.1038/s41467-020-15647-5CrossrefMedlineGoogle Scholar

• 42.

  Scharf GM, Kilian K, Cordero J, Wang Y, Grund A, Hofmann M, Froese N, Wang X, Kispert A, Kist R, et al.

  . Inactivation of Sox9 in fibroblasts reduces cardiac fibrosis and inflammation.JCI Insight2019; 4:e126721. doi: 10.1172/jci.insight.126721CrossrefGoogle Scholar

• 43.

  Liu G, Ma C, Yang H, Zhang PY

  . Transforming growth factor beta and its role in heart disease.Exp Ther Med2017; 13:2123–2128. doi: 10.3892/etm.2017.4246CrossrefMedlineGoogle Scholar

• 44.

  Dobaczewski M, Chen W, Frangogiannis NG

  . Transforming growth factor (TGF)-beta signaling in cardiac remodeling.J Mol Cell Cardiol2011; 51:600–606. doi: 10.1016/j.yjmcc.2010.10.033CrossrefMedlineGoogle Scholar

• 45.

  Frangogiannis NG, Ren G, Dewald O, Zymek P, Haudek S, Koerting A, Winkelmann K, Michael LH, Lawler J, Entman ML

  . Critical role of endogenous thrombospondin-1 in preventing expansion of healing myocardial infarcts.Circulation. 2005; 111:2935–2942. doi: 10.1161/CIRCULATIONAHA.104.510354LinkGoogle Scholar

• 46.

  Lal H, Ahmad F, Zhou J, Yu JE, Vagnozzi RJ, Guo Y, Yu D, Tsai EJ, Woodgett J, Gao E, et al.

  . Cardiac fibroblast glycogen synthase kinase-3β regulates ventricular remodeling and dysfunction in ischemic heart.Circulation. 2014; 130:419–430. doi: 10.1161/CIRCULATIONAHA.113.008364LinkGoogle Scholar

• 47.

  Bujak M, Ren G, Kweon HJ, Dobaczewski M, Reddy A, Taffet G, Wang XF, Frangogiannis NG

  . Essential role of Smad3 in infarct healing and in the pathogenesis of cardiac remodeling.Circulation. 2007; 116:2127–2138. doi: 10.1161/CIRCULATIONAHA.107.704197LinkGoogle Scholar

• 48.

  Campa CC, Silva RL, Margaria JP, Pirali T, Mattos MS, Kraemer LR, Reis DC, Grosa G, Copperi F, Dalmarco EM, et al.

  . Inhalation of the prodrug PI3K inhibitor CL27c improves lung function in asthma and fibrosis.Nat Commun. 2018; 9:5232. doi: 10.1038/s41467-018-07698-6CrossrefMedlineGoogle Scholar

• 49.

  Li J, Cao J, Li M, Yu Y, Yang Y, Xiao X, Wu Z, Wang L, Tu Y, Chen H

  . Collagen triple helix repeat containing-1 inhibits transforming growth factor-b1-induced collagen type I expression in keloid.Br J Dermatol2011; 164:1030–1036. doi: 10.1111/j.1365-2133.2011.10215.xCrossrefMedlineGoogle Scholar

• 50.

  Oka T, Xu J, Kaiser RA, Melendez J, Hambleton M, Sargent MA, Lorts A, Brunskill EW, Dorn GW, Conway SJ, et al.

  . Genetic manipulation of periostin expression reveals a role in cardiac hypertrophy and ventricular remodeling.Circ Res. 2007; 101:313–321. doi: 10.1161/CIRCRESAHA.107.149047LinkGoogle Scholar

• 51.

  Kaur H, Takefuji M, Ngai CY, Carvalho J, Bayer J, Wietelmann A, Poetsch A, Hoelper S, Conway SJ, Möllmann H, et al.

  . Targeted ablation of periostin-expressing activated fibroblasts prevents adverse cardiac remodeling in mice.Circ Res. 2016; 118:1906–1917. doi: 10.1161/CIRCRESAHA.116.308643LinkGoogle Scholar

• 52.

  Kim GH, Uriel N, Burkhoff D

  . Reverse remodelling and myocardial recovery in heart failure.Nat Rev Cardiol. 2018; 15:83–96. doi: 10.1038/nrcardio.2017.139CrossrefMedlineGoogle Scholar

• 53.

  Buch T, Heppner FL, Tertilt C, Heinen TJ, Kremer M, Wunderlich FT, Jung S, Waisman A

  . A Cre-inducible diphtheria toxin receptor mediates cell lineage ablation after toxin administration.Nat Methods. 2005; 2:419–426. doi: 10.1038/nmeth762CrossrefMedlineGoogle Scholar

• 54.

  Döring H

  . The isolated perfused heart according to Langendorff technique-function-application.Physiol Bohemoslov1990; 39:481–504MedlineGoogle Scholar

• 55.

  Jaitin DA, Kenigsberg E, Keren-Shaul H, Elefant N, Paul F, Zaretsky I, Mildner A, Cohen N, Jung S, Tanay A, et al.

  . Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types.Science. 2014; 343:776–779. doi: 10.1126/science.1247651CrossrefMedlineGoogle Scholar

• 56.

  Lavin Y, Kobayashi S, Leader A, Amir ED, Elefant N, Bigenwald C, Remark R, Sweeney R, Becker CD, Levine JH, et al.

  . Innate immune landscape in early lung adenocarcinoma by paired single-cell analyses.Cell. 2017; 169:750–765. doi: 10.1016/j.cell.2017.04.014CrossrefMedlineGoogle Scholar

• 57.

  Bagnoli JW, Ziegenhain C, Janjic A, Wange LE, Vieth B, Parekh S, Geuder J, Hellmann I, Enard W

  . Sensitive and powerful single-cell RNA sequencing using mcSCRB-seq.Nat Commun. 2018; 9:2937. doi: 10.1038/s41467-018-05347-6CrossrefMedlineGoogle Scholar

• 58.

  Lara-Astiaso D, Weiner A, Lorenzo-Vivas E, Zaretsky I, Jaitin DA, David E, Keren-Shaul H, Mildner A, Winter D, Jung S, et al.

  . Immunogenetics: chromatin state dynamics during blood formation.Science. 2014; 345:943–949. doi: 10.1126/science.1256271CrossrefMedlineGoogle Scholar

• 59.

  Corces MR, Buenrostro JD, Wu B, Greenside PG, Chan SM, Koenig JL, Snyder MP, Pritchard JK, Kundaje A, Greenleaf WJ, et al.

  . Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution.Nat Genet. 2016; 48:1193–1203. doi: 10.1038/ng.3646CrossrefMedlineGoogle Scholar

• 60.

  Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ

  . Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position.Nat Methods. 2013; 10:1213–1218. doi: 10.1038/nmeth.2688CrossrefMedlineGoogle Scholar

• 61.

  Stohn JP, Wang Q, Siviski ME, Kennedy K, Jin YR, Kacer D, DeMambro V, Liaw L, Vary CP, Rosen CJ

  . Cthrc1 controls adipose tissue formation, body composition, and physical activity.Obesity2015; 23:1633–1642. doi: 10.1002/oby.21144CrossrefMedlineGoogle Scholar

• 62.

  Ruiz-Villalba A, van Pelt-Verkuil E, Gunst QD, Ruijter JM, van den Hoff MJ

  . Amplification of nonspecific products in quantitative polymerase chain reactions (qPCR).Biomol Detect Quantif. 2017; 14:7–18. doi: 10.1016/j.bdq.2017.10.001CrossrefMedlineGoogle Scholar

• 63.

  Ruijter JM, Ramakers C, Hoogaars WM, Karlen Y, Bakker O, van den Hoff MJ, Moorman AF

  . Amplification efficiency: linking baseline and bias in the analysis of quantitative PCR data.Nucleic Acids Res. 2009; 37:e45. doi: 10.1093/nar/gkp045CrossrefMedlineGoogle Scholar

• 64.

  Ruijter JM, Ruiz Villalba A, Hellemans J, Untergasser A, van den Hoff MJ

  . Removal of between-run variation in a multi-plate qPCR experiment.Biomol Detect Quantif. 2015; 5:10–14. doi: 10.1016/j.bdq.2015.07.001CrossrefMedlineGoogle Scholar

• 65.

  Ruiz-Villalba A, Mattiotti A, Gunst QD, Cano-Ballesteros S, van den Hoff MJ, Ruijter JM

  . Reference genes for gene expression studies in the mouse heart.Sci Rep. 2017; 7:24. doi: 10.1038/s41598-017-00043-9CrossrefMedlineGoogle Scholar

• 66.

  Nygard AB, Jørgensen CB, Cirera S, Fredholm M

  . Selection of reference genes for gene expression studies in pig tissues using SYBR green qPCR.BMC Mol Biol. 2007; 8:67. doi: 10.1186/1471-2199-8-67CrossrefMedlineGoogle Scholar

• 67.

  Villani A-C, Satija R, Reynolds G, Sarkizova S, Shekhar K, Fletcher J, Griesbeck M, Butler A, Zheng S, Lazo S

  . Single-cell RNA-seq reveals new types of human blood dendritic cells, monocytes, and progenitors.Science2017; 356:eaah4573. doi: 10.1126/science.aah4573CrossrefMedlineGoogle Scholar

• 68.

  Butler A, Hoffman P, Smibert P, Papalexi E, Satija R

  . Integrating single-cell transcriptomic data across different conditions, technologies, and species.Nat Biotechnol. 2018; 36:411–420. doi: 10.1038/nbt.4096CrossrefMedlineGoogle Scholar

• 69.

  Li W, Cerise JE, Yang Y, Han H

  . Application of t-SNE to human genetic data.J Bioinform Comput Biol. 2017; 15:1750017. doi: 10.1142/S0219720017500172CrossrefMedlineGoogle Scholar

• 70.

  Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR

  . STAR: ultrafast universal RNA-seq aligner.Bioinformatics. 2013; 29:15–21. doi: 10.1093/bioinformatics/bts635CrossrefMedlineGoogle Scholar

• 71. R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria; 2018. https://www./Google Scholar

• 72.

  Liao Y, Smyth GK, Shi W

  . The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads.Nucleic Acids Res2019; 47:e47. doi: 10.1093/nar/gkz114CrossrefMedlineGoogle Scholar

• 73.

  Robinson MD, McCarthy DJ, Smyth GK

  . edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.Bioinformatics. 2010; 26:139–140. doi: 10.1093/bioinformatics/btp616CrossrefMedlineGoogle Scholar

• 74.

  Love MI, Huber W, Anders S

  . Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2Genome Biol2014; 15:550. doi: 10.1186/s13059-014-0550-8CrossrefMedlineGoogle Scholar

• 75.

  Yu G, Wang LG, Han Y, He QY

  . clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS. 2012; 16:284–287. doi: 10.1089/omi.2011.0118CrossrefMedlineGoogle Scholar

• 76.

  Fabregat A, Jupe S, Matthews L, Sidiropoulos K, Gillespie M, Garapati P, Haw R, Jassal B, Korninger F, May B, et al.

  . The reactome pathway knowledgebase.Nucleic Acids Res. 2018; 46:D649–D655. doi: 10.1093/nar/gkx1132CrossrefMedlineGoogle Scholar

• 77.

  Ferreira JA

  . The Benjamini-Hochberg method in the case of discrete test statistics.Int J Biostat. 2007; 3:Article 11. doi: 10.2202/1557-4679.1065CrossrefMedlineGoogle Scholar

• 78.

  Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, et al.

  . Enrichr: a comprehensive gene set enrichment analysis web server 2016 update.Nucleic Acids Res. 2016; 44:W90–W97. doi: 10.1093/nar/gkw377CrossrefMedlineGoogle Scholar

• 79.

  Huang R, Grishagin I, Wang Y, Zhao T, Greene J, Obenauer JC, Ngan D, Nguyen D-T, Guha R, Jadhav A

  . The NCATS BioPlanet: an integrated platform for exploring the universe of cellular signaling pathways for toxicology, systems biology, and chemical genomics.Front Pharmacol2019; 10:445. doi: 10.3389/fphar.2019.00445CrossrefMedlineGoogle Scholar

• 80.

  Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, Hao Y, Stoeckius M, Smibert P, Satija R

  . Comprehensive integration of single-cell data.Cell. 2019; 177:1888–1902. doi: 10.1016/j.cell.2019.05.031CrossrefMedlineGoogle Scholar

• 81.

  Rouillard AD, Gundersen GW, Fernandez NF, Wang Z, Monteiro CD, McDermott MG, Ma’ayan A

  . The harmonizome: a collection of processed datasets gathered to serve and mine knowledge about genes and proteins.Database2016; 2016:baw100. doi: 10.1093/database/baw100CrossrefMedlineGoogle Scholar

• 82.

  Han H, Cho JW, Lee S, Yun A, Kim H, Bae D, Yang S, Kim CY, Lee M, Kim E, et al.

  . TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions.Nucleic Acids Res. 2018; 46:D380–D386. doi: 10.1093/nar/gkx1013CrossrefMedlineGoogle Scholar

• 83.

  Perfetto L, Briganti L, Calderone A, Cerquone Perpetuini A, Iannuccelli M, Langone F, Licata L, Marinkovic M, Mattioni A, Pavlidou T, et al.

  . SIGNOR: a database of causal relationships between biological entities.Nucleic Acids Res. 2016; 44:D548–D554. doi: 10.1093/nar/gkv1048CrossrefMedlineGoogle Scholar

• 84. The UniProt Consortium. UniProt: the universal protein knowledgebase.Nucleic Acids Res. 2017; 45:D158–D169. doi: 10.1093/nar/gkw1099CrossrefMedlineGoogle Scholar

• 85.

  Langmead B, Trapnell C, Pop M, Salzberg SL

  . Ultrafast and memory-efficient alignment of short DNA sequences to the human genome.Genome Biol. 2009; 10:R25. doi: 10.1186/gb-2009-10-3-r25CrossrefMedlineGoogle Scholar

• 86.

  Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al.

  . Model-based analysis of ChIP-Seq (MACS).Genome Biol. 2008; 9:R137. doi: 10.1186/gb-2008-9-9-r137CrossrefMedlineGoogle Scholar

• 87.

  Yu G, Wang LG, He QY

  . ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization.Bioinformatics. 2015; 31:2382–2383. doi: 10.1093/bioinformatics/btv145CrossrefMedlineGoogle Scholar

• 88.

  Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK

  . Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities.Mol Cell. 2010; 38:576–589. doi: 10.1016/j.molcel.2010.05.004CrossrefMedlineGoogle Scholar

• 89.

  Scrucca L, Fop M, Murphy TB, Raftery AE

  . mclust 5: Clustering, classification and density estimation using gaussian finite mixture models.R J. 2016; 8:289–317.CrossrefMedlineGoogle Scholar