Advertisement
Journal of Investigative Dermatology Home

Single-Cell RNA-Sequencing Reveals Lineage-Specific Regulatory Changes of Fibroblasts and Vascular Endothelial Cells in Keloids

  • Author Footnotes
    3 These authors contributed equally.
    Xuanyu Liu
    Footnotes
    3 These authors contributed equally.
    Affiliations
    State Key Laboratory of Cardiovascular Disease, Beijing Key Laboratory for Molecular Diagnostics of Cardiovascular Diseases, Center of Laboratory Medicine, Fuwai Hospital, National Center for Cardiovascular Diseases, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
    Search for articles by this author
  • Author Footnotes
    3 These authors contributed equally.
    Wen Chen
    Footnotes
    3 These authors contributed equally.
    Affiliations
    State Key Laboratory of Cardiovascular Disease, Beijing Key Laboratory for Molecular Diagnostics of Cardiovascular Diseases, Center of Laboratory Medicine, Fuwai Hospital, National Center for Cardiovascular Diseases, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
    Search for articles by this author
  • Author Footnotes
    3 These authors contributed equally.
    Qingyi Zeng
    Footnotes
    3 These authors contributed equally.
    Affiliations
    State Key Laboratory of Cardiovascular Disease, Beijing Key Laboratory for Molecular Diagnostics of Cardiovascular Diseases, Center of Laboratory Medicine, Fuwai Hospital, National Center for Cardiovascular Diseases, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
    Search for articles by this author
  • Baihui Ma
    Affiliations
    State Key Laboratory of Cardiovascular Disease, Beijing Key Laboratory for Molecular Diagnostics of Cardiovascular Diseases, Center of Laboratory Medicine, Fuwai Hospital, National Center for Cardiovascular Diseases, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
    Search for articles by this author
  • Zhujun Li
    Affiliations
    Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing, China
    Search for articles by this author
  • Tian Meng
    Affiliations
    Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing, China
    Search for articles by this author
  • Jie Chen
    Affiliations
    Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing, China
    Search for articles by this author
  • Nanze Yu
    Affiliations
    Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing, China
    Search for articles by this author
  • Zhou Zhou
    Affiliations
    State Key Laboratory of Cardiovascular Disease, Beijing Key Laboratory for Molecular Diagnostics of Cardiovascular Diseases, Center of Laboratory Medicine, Fuwai Hospital, National Center for Cardiovascular Diseases, Chinese Academy of Medical Sciences and Peking Union Medical College, Beijing, China
    Search for articles by this author
  • Xiao Long
    Correspondence
    Correspondence: Xiao Long, Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing 100730, China.
    Affiliations
    Division of Plastic Surgery, Peking Union Medical College Hospital, Beijing, China
    Search for articles by this author
  • Author Footnotes
    3 These authors contributed equally.
Open AccessPublished:July 06, 2021DOI:https://doi.org/10.1016/j.jid.2021.06.010
      Keloids are a benign dermal fibrotic disorder with features similar to malignant tumors. keloids remain a therapeutic challenge and lack medical therapies, which is partially due to the incomplete understanding of the pathogenesis mechanism. We performed single-cell RNA-sequencing of 28,064 cells from keloid skin tissue and adjacent relatively normal tissue. Unbiased clustering revealed substantial cellular heterogeneity of keloid tissue, which included 21 clusters assigned to 11 cell lineages. We observed significant expansion of fibroblast and vascular endothelial cell subpopulations in keloids, reflecting their strong association with keloid pathogenesis. Comparative analyses were performed to identify the dysregulated pathways, regulators and ligand-receptor interactions in keloid fibroblasts and vascular endothelial cells. Our results highlight the roles of TGFβ and Eph-ephrin signaling pathways in both the aberrant fibrogenesis and angiogenesis of keloids. Critical regulators probably involved in the fibrogenesis of keloid fibroblasts, such as TWIST1, FOXO3 and SMAD3, were identified. TWIST1 inhibitor harmine could significantly suppress the fibrogenesis of keloid fibroblasts. In addition, tumor-related pathways were activated in keloid fibroblasts and vascular endothelial cells, which may be responsible for the malignant features of keloids. Our study put insights into the pathogenesis of keloids and provides potential targets for medical therapies.

      Graphical abstract

      Abbreviations:

      ECM (extracellular matrix), RNA-seq (RNA-sequencing), TF (transcription factor), VEC (vascular endothelial cell)

      Introduction

      Keloids are a dermal fibrotic disorder arising following an aberrant wound-healing response (
      • Glass 2nd, D.A.
      Current understanding of the genetic causes of keloid formation.
      ). Histologically, keloid scars are characterized by excessive extracellular matrix (ECM) deposition and a rich vasculature (
      • Ashcroft K.J.
      • Syed F.
      • Bayat A.
      Site-specific keloid fibroblasts alter the behaviour of normal skin and normal scar fibroblasts through paracrine signalling.
      ). Although keloids are classified as a benign dermal growth, they exhibit biological features similar to malignant tumors, such as hyperproliferation, apoptosis resistance and invasion (
      • Mari W.
      • Alsabri S.G.
      • Tabal N.
      • Younes S.
      • Sherif A.
      • Simman R.
      Novel insights on understanding of keloid scar: article review.
      ). Keloid scars grow beyond the boundaries of the original wound, causing pain, pruritus and contracture, which leads to a serious physical and psychological burden for patients (
      • Gauglitz G.G.
      • Korting H.C.
      • Pavicic T.
      • Ruzicka T.
      • Jeschke M.G.
      Hypertrophic scarring and keloids: pathomechanisms and current and emerging treatment strategies.
      ). Although a wide range of therapies are currently being used, keloids remain a therapeutic challenge in terms of high recurrence rate and lack of medical therapies (
      • Mari W.
      • Alsabri S.G.
      • Tabal N.
      • Younes S.
      • Sherif A.
      • Simman R.
      Novel insights on understanding of keloid scar: article review.
      ). This is partially due to the incomplete understanding of keloid pathogenesis, although both environmental and genetic factors have been implicated (
      • Nakashima M.
      • Chung S.
      • Takahashi A.
      • Kamatani N.
      • Kawaguchi T.
      • Tsunoda T.
      • et al.
      A genome-wide association study identifies four susceptibility loci for keloid in the Japanese population.
      ;
      • Shih B.
      • Garside E.
      • McGrouther D.A.
      • Bayat A.
      Molecular dissection of abnormal wound healing processes resulting in keloid disease.
      ). A thorough understanding of the cellular and molecular mechanisms of keloid pathogenesis would facilitate the development of medical therapies for this disease.
      Recent technical advances in single-cell RNA-sequencing (RNA-seq) have enabled the transcriptomes of tens of thousands of cells to be assayed at a single-cell resolution (
      • Zheng G.X.
      • Terry J.M.
      • Belgrader P.
      • Ryvkin P.
      • Bent Z.W.
      • Wilson R.
      • et al.
      Massively parallel digital transcriptional profiling of single cells.
      ). Compared with the average expression of genes from a mixed cell population obtained via bulk RNA-seq, large-scale single-cell RNA-seq allows unbiased assessment of cellular heterogeneity and regulatory network at an unprecedented scale and resolution (
      • Kulkarni A.
      • Anderson A.G.
      • Merullo D.P.
      • Konopka G.
      Beyond bulk: a review of single cell transcriptomics methodologies and applications.
      ). Single-cell RNA-seq is therefore emerging as a powerful tool for understanding the cellular and molecular mechanisms of the pathogenesis of a variety of diseases, including fibrotic disorders such as pulmonary fibrosis (
      • Reyfman P.A.
      • Walter J.M.
      • Joshi N.
      • Anekalla K.R.
      • McQuattie-Pimentel A.C.
      • Chiu S.
      • et al.
      Single-cell transcriptomic analysis of human lung provides insights into the pathobiology of pulmonary fibrosis.
      ) and lupus nephritis (
      • Der E.
      • Suryawanshi H.
      • Morozov P.
      • Kustagi M.
      • Goilav B.
      • Ranabathou S.
      • et al.
      Tubular cell and keratinocyte single-cell transcriptomics applied to lupus nephritis reveal type I IFN and fibrosis relevant pathways [published correction appears in Nat Immunol 2019;20:1556].
      ). Previous efforts have been made to examine the transcriptomic alterations in keloid tissue using bulk RNA-seq or microarray analysis (
      • Liang X.
      • Ma L.
      • Long X.
      • Wang X.
      LncRNA expression profiles and validation in keloid and normal skin tissue.
      ;
      • Onoufriadis A.
      • Hsu C.K.
      • Ainali C.
      • Ung C.Y.
      • Rashidghamat E.
      • Yang H.S.
      • et al.
      Time series integrative analysis of RNA sequencing and microRNA expression data reveals key biologic wound healing pathways in keloid-prone individuals.
      ;
      • Wang J.
      • Wu H.
      • Xiao Z.
      • Dong X.
      Expression profiles of lncRNAs and circRNAs in keloid.
      ). However, to our knowledge, the cellular heterogeneity and regulatory changes in keloids have not yet been systematically investigated at a single-cell resolution.
      In this study, we performed single-cell RNA-seq of 28,064 cells from keloid skin tissue and adjacent normal tissue. Comparative analyses were performed to identify the dysregulated pathways, regulators and ligand-receptor interactions in keloid fibroblasts and vascular endothelial cells (VECs), the two important cell lineages in keloid pathogenesis and medical interventions. Our study put insights into the pathogenesis of keloids and provided potential targets for medical therapies. Our dataset also constitutes a valuable resource for further investigations of the mechanism of keloid pathogenesis.

      Results

       Single-cell RNA-seq reveals the cellular diversity and heterogeneity of keloid skin tissue

      To dissect cellular heterogeneity and explore the regulatory changes in keloid skin tissue, we sampled keloid lesional skin tissue (i.e., CASE) and matched normal tissue adjacent to the keloid scar (i.e., CTRL) from four patients (Supplementary Table S1). The eight samples were dissociated into single cells and subjected to single-cell RNA-seq (Figure 1a). After quality control, the transcriptomes of 28,064 cells (CASE: 12,425; CTRL: 15,639) were obtained. Unbiased clustering revealed 21 cell clusters (Figure 1b). Based on hierarchical clustering (Figure 1c) and established lineage-specific marker genes (Figure 1d), we assigned these clusters to 11 cell lineages. The keratinocyte lineage (marked by KRT5 and KRT14) (
      • Joost S.
      • Zeisel A.
      • Jacob T.
      • Sun X.
      • La Manno G.
      • Lönnerberg P.
      • et al.
      Single-cell transcriptomics reveals that differentiation and spatial signatures shape epidermal and hair follicle heterogeneity.
      ) accounted for the largest proportion (41.7%) of cells. VECs (marked by PECAM1 and CDH5) (
      • Kalucka J.
      • de Rooij L.P.M.H.
      • Goveia J.
      • Rohlenova K.
      • Dumas S.J.
      • Meta E.
      • et al.
      Single-cell transcriptome atlas of murine endothelial cells.
      ) accounted for the second-largest proportion (26.3%) of cells. Lymphatic endothelial cells (c16) differ from VECs in the expression of lineage-specific markers such as LYVE1 and PROX1 (
      • Johnson N.C.
      • Dillard M.E.
      • Baluk P.
      • McDonald D.M.
      • Harvey N.L.
      • Frase S.L.
      • et al.
      Lymphatic endothelial cell identity is reversible and its maintenance requires Prox1 activity.
      ). The fibroblast lineage (marked by PDGFRA and DCN) (
      • Guerrero-Juarez C.F.
      • Dedhia P.H.
      • Jin S.
      • Ruiz-Vega R.
      • Ma D.
      • Liu Y.
      • et al.
      Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds.
      ) included two clusters c3 and c9. Also, we found other typical skin lineages including sweat gland cells (marked by AQP5 and MUCL1) (
      • He H.
      • Suryawanshi H.
      • Morozov P.
      • Gay-Mimbrera J.
      • Del Duca E.
      • Kim H.J.
      • et al.
      Single-cell transcriptome analysis of human skin identifies novel fibroblast subpopulation and enrichment of immune subsets in atopic dermatitis.
      ), smooth muscle cells (marked by MYH11 and CNN1) (
      • Liu X.
      • Chen W.
      • Li W.
      • Li Y.
      • Priest J.R.
      • Zhou B.
      • et al.
      Single-cell RNA-seq of the developing cardiac outflow tract reveals convergent development of the vascular smooth muscle cells.
      ), mural cells (marked by PDGFRB and RGS5) (
      • Holm A.
      • Heumann T.
      • Augustin H.G.
      Microvascular mural cell organotypic heterogeneity and functional plasticity.
      ), neural cells (marked by NRXN1 and SCN7A) (
      • He H.
      • Suryawanshi H.
      • Morozov P.
      • Gay-Mimbrera J.
      • Del Duca E.
      • Kim H.J.
      • et al.
      Single-cell transcriptome analysis of human skin identifies novel fibroblast subpopulation and enrichment of immune subsets in atopic dermatitis.
      ), melanocytes (marked by MLANA and DCT) (
      • Miller A.J.
      • Du J.
      • Rowan S.
      • Hershey C.L.
      • Widlund H.R.
      • Fisher D.E.
      Transcriptional regulation of the melanoma prognostic marker melastatin (TRPM1) by MITF in melanocytes and melanoma.
      ) and Schwann cells (marked by GLDN and TJP1) (
      • He H.
      • Suryawanshi H.
      • Morozov P.
      • Gay-Mimbrera J.
      • Del Duca E.
      • Kim H.J.
      • et al.
      Single-cell transcriptome analysis of human skin identifies novel fibroblast subpopulation and enrichment of immune subsets in atopic dermatitis.
      ). These clusters showed distinct molecular signatures (Figure 1e and Supplementary Table S2), reflecting the cellular diversity and heterogeneity of keloid skin tissues.
      Figure thumbnail gr1
      Figure 1Single-cell RNA-seq reveals the cellular diversity and heterogeneity of keloid skin tissue. (a) Schematic representation of the experimental procedure. CASE and CTRL skin tissues were harvested separately at the time of surgery (n = 4). (b) Unbiased clustering of 28,064 cells reveals 21 cellular clusters. Clusters are distinguished by different colors. The number in the parenthesis indicates the cell count. (c) Hierarchical clustering of the clusters based on the average expression of the 2,000 most variable genes. (d) Expression of the established marker genes specific for each lineage in each cluster. (e) Representative molecular signatures for each cell cluster. The area of the circles indicates the proportion of cells expressing the gene, and the color intensity reflects the expression intensity. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FB, fibroblast; KTR, keratinocyte; lEndo, lymphatic endothelial cell; LEU, Leukocyte; MLA, melanocyte; RNA-seq, RNA-sequencing; SGC, sweat gland cell; SMC, smooth muscle cell; UMAP, uniform manifold approximation and projection; vEndo, vascular endothelial cell.

       Differential proportion analysis reveals significant expansion of fibroblast and VEC subpopulations in keloids

      We next attempted to identify keloid-associated cell lineages or clusters that were significantly expanded or contracted in keloids. The visualization of cellular density and distribution revealed dramatic changes in the relative proportions of multiple cell lineages (Figure 2a and b and Supplementary Figure S1). For example, lineage expansion of VECs and fibroblasts was observed in keloids, while contraction was observed for keratinocytes and leukocytes. We further performed statistical tests, taking individual patients into account (Supplementary Figure S2 and Supplementary Table S3). Only the change in VECs reached statistical significance (paired t-test P-value = 0.01). Similarly, we performed statistical tests at the cluster level (Figure 2c). Three VEC clusters, c4, c5 and c18, were significantly expanded in CASE. Notably, one fibroblast cluster, c9, was also significantly expanded. Overall, differential proportion analysis revealed significant expansion of fibroblast and VEC subpopulations in keloids compared with normal skin tissue, reflecting their strong association with keloid pathogenesis. This study thus focused on these two cell lineages.
      Figure thumbnail gr2
      Figure 2Differential proportion analysis reveals significant expansion of fibroblast and vascular endothelial cell subpopulations in keloids compared with relatively normal skin tissue. (a) Visualization of cellular density reveals dramatic changes in the proportions of multiple cell lineages in CASE versus CTRL. (b) Distribution of CASE and CTRL cells in the UMAP embeddings. In a and b, equal number of cells were sampled in CASE (n = 12,425) and CTRL (n = 12,425). (c) Vascular endothelial clusters c4, c5 and c18 and fibroblast cluster c9 are significantly expanded in CASE versus CTRL. The relative proportion of each cluster was calculated in each sample. A significance threshold of P-value < 0.05 of paired t-test was used. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; KTR, keratinocyte; lEndo, lymphatic endothelial cell; LEU, Leukocyte; MLA, melanocyte; SGC, sweat gland cell; UMAP, uniform manifold approximation and projection.

       Gene set enrichment analysis reveals dysregulated pathways in keloid fibroblasts

      We next identified fibroblast-specific differentially regulated pathways in keloid versus normal skin tissue (Figure 3a and Supplementary Table S4) through gene set enrichment analysis. Consistent with the excessive ECM deposition observed in keloids, ECM-related pathways, such as ECM organization and collagen formation, were significantly upregulated (Figure 3a; gene set enrichment analysis false discovery rate q-value < 0.05). Notably, several signaling pathways were activated in keloid fibroblasts (Figure 3b). The role of TGFβ and canonical WNT signaling pathways in keloids or other fibrotic disorders has been well established (
      • Piersma B.
      • Bank R.A.
      • Boersema M.
      Signaling in fibrosis: TGF-β, WNT, and YAP/TAZ converge.
      ). In line with this, we observed the activation of these two master pathways (Figure 3b). PTEN, a known tumor suppressor, is a major growth signaling inhibitor that controls cell growth, survival and proliferation (
      • Chalhoub N.
      • Baker S.J.
      PTEN and the PI3-kinase pathway in cancer.
      ). Intriguingly, we found that a set of genes involved in the negative regulation of PTEN gene transcription were upregulated in keloid fibroblasts (Figure 3b), which agrees with the malignant features of keloid fibroblasts such as excessive proliferation and resistance to apoptosis (
      • Lim C.P.
      • Phan T.T.
      • Lim I.J.
      • Cao X.
      Stat3 contributes to keloid pathogenesis via promoting collagen production, cell proliferation and migration.
      ). Notably, the Eph-ephrin signaling pathway, which has recently been recognized to be associated with cardiac fibrosis (
      • Su S.A.
      • Yang D.
      • Wu Y.
      • Xie Y.
      • Zhu W.
      • Cai Z.
      • et al.
      EphrinB2 regulates cardiac fibrosis through modulating the interaction of Stat3 and TGF-β/Smad3 signaling.
      ), was significantly activated in keloid fibroblasts (Figure 3b). Multiple genes encoding the ligands and receptors of Eph-ephrin signaling, such as EFNB1, EFNB2, EPHB2, EPHB3 and EPHA3 were upregulated in keloid fibroblasts.
      Figure thumbnail gr3
      Figure 3Gene set enrichment analysis reveals dysregulated pathways in keloid fibroblasts. (a) Network view of differentially regulated REACTOME pathways in keloid fibroblasts. The size of the dot reflects the size of the gene set. The dots in red denote upregulated pathways, and dots in blue represent downregulated pathways. A significance threshold of an FDR q-value of 0.05 was used. (b) Enrichment plots (upper panel) and leading-edge gene expression heatmaps (the top 20 genes; lower panel) for representative signaling pathways upregulated in keloid fibroblasts. NES used to compare analysis results across gene sets. The vertical lines in the enrichment plot show where the members of the gene set appear in the ranked list of genes. Leading-edge genes: the subset of genes in the gene set that contribute most to the enrichment. The average expression across cells in each group is shown in the heatmap. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FDR, false discovery rate; NES, normalized enrichment score.

       Cellular heterogeneity of keloid fibroblasts

      We next dissected the cellular heterogeneity and characterized the subpopulations of fibroblasts. Fibroblast clusters c3 and c9 showed distinct expression profiles (Figure 4a). Compared with c3, c9 displayed more obvious phenotypes of myofibroblasts. For example, significantly higher expression of contractile genes (e.g., ACTA2, TAGLN and MYL9), collagen and elastin genes (e.g., COL1A2, COL3A1 and ELN), as well as myofibroblast markers (e.g., FN1 and CTHRC1), was observed (Figure 4b; adjusted P-value < 0.01). This finding agrees with the significant expansion of c9 in keloids (Figure 2c), reflecting the elevated myofibrogenesis. To further dissect the cellular heterogeneity, we performed secondary clustering and identified five subpopulations (Figure 4c), which exhibited distinct molecular signatures (Figure 4d and Supplementary Table S5). These subpopulations could be marked by specific markers: s0-PDGFRA+ POSTNhigh, s1-PDGFRA+ ABCA8high, s2-PDGFRA+ RELhigh, s3-PDGFRA+ RGS2high and s4-PDGFRA+ WISP2high (Figure 4e). Correlation analysis was performed to infer the relationships between them (Figure 4f). For example, the closest subpopulation to s0 was s4. Functional enrichment analysis revealed that the signature genes of s0 and s4 were both enriched with ECM-related terms, thus representing two myofibroblast subpopulations (Figure 4g). However, they differed in the components of the ECM: s0 exhibited a higher level of collagens, and s4 exhibited a higher level of elastin (Figure 4d). The signature genes of s1 were enriched with terms related to immune responses, such as the regulation of complement activation and lymphocyte chemotaxis. The signature genes of s2 were enriched with terms related to leukocyte chemotaxis, carbohydrate metabolism and stress responses. The signature genes of s3 were enriched with terms related to stress responses and differentiation. Immunofluorescence staining results confirmed the presence of the five subpopulations in keloid tissue (Figure 4h). Supplementary Figure S3 shows low-magnification overviews of the staining.
      Figure thumbnail gr4
      Figure 4Cellular heterogeneity of keloid fibroblasts revealed by single-cell transcriptomic data. (a) Heatmap showing distinct expression profiles between the fibroblast clusters c3 and c9. (b) Compared with c3, c9 displays more obvious phenotypes of myofibroblasts. ∗∗: adjusted P-value < 0.01. (c) Secondary clustering of fibroblasts further identifies five subpopulations. The upper small panel shows the expression distribution of ACTA2. The lower small panel shows the distribution of CASE and CTRL cells. (d) The five fibroblast subpopulations display distinct expression profiles. (e) Molecular signatures for the five fibroblast subpopulations. (f) Correlations among the five subpopulations. (g) Functional enrichment for each of the five subpopulations. The significance threshold was set to an adjusted P-value < 0.05. (h) The presence of the five subpopulations of keloid fibroblasts was confirmed in keloid tissue by immunofluorescence staining. Bar = 5 μm. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; STAT, signal transducer and activator of transcription.

       Pseudotemporal ordering of fibroblasts reveals a branched trajectory with a significant shift toward the myofibroblast phenotype in keloids

      To further explore the relationships among the fibroblast subpopulations and study the regulatory dynamics during the fibroblast-to-myofibroblast phenotypic transition, we performed pseudotemporal ordering of all fibroblasts with Monocle2 (
      • Qiu X.
      • Hill A.
      • Packer J.
      • Lin D.
      • Ma Y.A.
      • Trapnell C.
      Single-cell mRNA quantification and differential analysis with Census.
      ). This analysis revealed a branched trajectory with two major branches, that is, cell fate 1 and cell fate 2 (Figure 5a). Subpopulations s0 and s4 constituted the most of the cell fate1 branch, which thus represents the cellular states of the myofibroblast phenotype (Figure 5b). Subpopulation s1 and s3 constituted the most of the cell fate 2 branch, which thus represents special states of cells associated with the immune response, stress response and differentiation. Subpopulation s2 accounted for the largest proportion of the prebranch, which represents the initial states of fibroblasts. Notably, compared with the normal control, the trajectory in keloids displayed a significant shift toward the myofibroblast phenotype (fate 1 cells accounted for 28% of the cells in CTRL versus 40% of those in CASE; Figure 5c). Through branched expression analysis modeling tests, we obtained the expression dynamics of 1,204 branch-dependent genes during the cellular state transition from the prebranch to cell fate 1 and cell fate 2 (Figure 5d and Supplementary Table S6; q-value < 1E–04). The hierarchical clustering of these genes revealed five gene modules, and the representative transcription factors (TFs) for each gene module are shown in Figure 5d. Among these modules, most of the genes in the module I exhibited high expression levels in the prebranch cells. Notably, TWIST1 is a representative TF of the module I, which has recently been recognized as an important profibrotic regulator in a variety of organ fibrotic disorders, including skin fibrosis (
      • Ning X.
      • Zhang K.
      • Wu Q.
      • Liu M.
      • Sun S.
      Emerging role of Twist1 in fibrotic diseases.
      ).
      Figure thumbnail gr5
      Figure 5Pseudo-temporal ordering and gene regulatory network analysis of keloid fibroblasts. (a) Pseudo-temporal ordering of keloid fibroblasts reveals a branched trajectory. (b) Distribution of the five subpopulations on each of the branches. (c) Significant shift toward a myofibroblast phenotype (cell fate 1) in CASE (right panel) versus CTRL (left panel). (d) Hierarchical clustering of the branch-dependent genes reveals five gene modules. The significant threshold was set to a q-value of the branched expression analysis modeling test < 1E-04. Representative TFs are shown. (e) Comparative analysis of the gene regulatory networks of fibroblasts between CASE (right panel) and CTRL (left panel) reveals dysregulated genes in keloid fibroblasts. The node size reflects degree centrality. The TFs (in black) and receptors (in orange) in the top 100 genes ranked by the delta degree are shown. (f) The nine TFs exhibiting both branch-dependence and dysregulation in the keloid network are all from gene module I, most of which are highly expressed on the prebranch. (g) Split violin plot showing the expression of the nine TFs in CASE and CTRL. BEAM, branched expression analysis modeling; CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; TF, transcription factor.

       Comparative analysis of the gene regulatory networks reveals dysregulated genes in keloid fibroblasts

      To further prioritize the TFs identified above, we built gene regulatory networks from single-cell data using a method implemented in bigScale2 (
      • Iacono G.
      • Massoni-Badosa R.
      • Heyn H.
      Single-cell transcriptomics unveils gene regulatory network plasticity.
      ), which allowed us to quantify the biological importance of genes and find the key regulators under diseased conditions. Figure 5e shows the regulatory networks constructed for fibroblasts from normal and keloid skin tissue. Comparative analysis between the CASE and CTRL networks revealed a group of genes that were greatly altered in terms of degree centrality (the number of edges afferent to a given node; Supplementary Table S7). We focused on the TFs and signal receptors of the top 100 genes that were dysregulated in keloids, which are indicated in Figure 5e and f. Notably, among the TF genes, TWIST1 was ranked at the top according to degree centrality, suggesting a key role of this regulator in keloid pathogenesis. EPHB2 was ranked at the top among the signal receptor genes, consolidating our view that the Eph-ephrin signaling pathway may play critical roles in keloid pathogenesis. We next tried to identify TFs that were both branch-dependent (Figure 5d) and dysregulated in the keloid network (Figure 5e), which would be more reliable given the independent analysis approaches. We ultimately identified nine TFs including TWIST1, YBX3, NFIL3, JARID2, NCOA7, FOXO3, NFKBIA, ETS2 and SMAD3 (Figure 5f). Among these TFs, SMAD3 is a known regulator of keloid through which TGFβ signaling exerts its profibrotic effects (
      • Darby I.A.
      • Laverdet B.
      • Bonté F.
      • Desmoulière A.
      Fibroblasts and myofibroblasts in wound healing.
      ). TWIST1 and FOXO3 have been implicated in fibrogenesis and organ fibrosis (
      • Al-Tamari H.M.
      • Dabral S.
      • Schmall A.
      • Sarvari P.
      • Ruppert C.
      • Paik J.
      • et al.
      FoxO3 an important player in fibrogenesis and therapeutic target for idiopathic pulmonary fibrosis.
      ;
      • Ning X.
      • Zhang K.
      • Wu Q.
      • Liu M.
      • Sun S.
      Emerging role of Twist1 in fibrotic diseases.
      ). Intriguingly, we found that all nine candidate regulators were from the gene module I, most of which were highly expressed in the prebranch (Figure 5f). We next examined the expression of these regulators in fibroblasts from keloid and normal tissues (Figure 5g). The expression modes of these genes agreed well with those reported in other fibrotic disorders. For example, TWIST1 and SMAD3 were upregulated, and FOXO3 was downregulated in keloid versus normal fibroblasts.

       TWIST1 inhibitor harmine significantly suppresses the fibrogenesis of keloid fibroblasts

      To experimentally confirm the role of TWIST1 in keloid pathogenesis, we investigated the effects of TWIST1 inhibitor harmine (
      • Yochum Z.A.
      • Cades J.
      • Mazzacurati L.
      • Neumann N.M.
      • Khetarpal S.K.
      • Chatterjee S.
      • et al.
      A first-in-class twist1 inhibitor with activity in oncogene-driven lung cancer.
      ) on the fibrogenesis of keloid fibroblasts in vitro. Keloid and normal fibroblasts were treated with harmine at concentrations of 0 (DMSO), 10 and 20 μM for 48 hours (Figure 6a and b) separately. Differentially expression analysis of bulk RNA-seq data between 10 μM harmine and DMSO treatments revealed hundreds of genes whose expression was significantly altered in keloid or normal fibroblasts (Figure 6c and Supplementary Tables S8 and S9). Strikingly, ECM-related pathways, such as ECM organization, were significantly downregulated (q-value < 0.05, log2 fold change < -1) following 10 μM harmine treatment in keloid fibroblasts (Figure 6d), indicating that the fibrogenesis was significantly suppressed. In contrast, the downregulated genes in normal fibroblasts following harmine treatment were not significantly enriched for ECM-related pathways (Figure 6e). The mRNA expression of TWIST1 was significantly suppressed (q-value < 0.05) in both keloid and normal fibroblasts following 10 μM harmine treatment (Figure 6f and g). The expression of genes encoding established fibrotic factors, such as COL1A1, COL3A1 and FN1, were significantly suppressed in both keloid and normal fibroblasts. The expression of key genes in TGFβ signaling, such as TGFBR2 and TGFB1, was also downregulated in both keloid and normal fibroblasts. However, the fold changes of these downregulated genes in normal fibroblasts were lower compared with that in keloid fibroblasts. The expression of these genes could be confirmed by real-time PCR (Supplementary Figure S4 and Supplementary Table S10). Furthermore, immunoblot assays showed that the protein levels of TWIST1, Collagen-1, Collagen-3, SMAD2/3 were significantly decreased following harmine treatment in keloid fibroblasts (Figure 6h, P < 0.05). A reduced level of phosphorylated SMAD2/3 (pSMAD2/3) was observed in keloid fibroblasts, despite that it did not reach statistical significance. Compared with keloid fibroblasts, harmine treatment had little or moderate effects on normal fibroblasts at the protein level (Figure 6i), probably due to the relatively low expression of TWIST1 in normal fibroblasts. Together, our results supported the role of TWIST1 in the fibrogenesis of keloid fibroblasts, and harmine may serve as a targeted drug for keloid treatment.
      Figure thumbnail gr6
      Figure 6TWIST1 inhibitor harmine significantly suppresses the fibrogenesis of keloid fibroblasts. (a) Schematic representation of the experimental procedure. Keloid or normal fibroblasts were isolated and treated with TWIST1 inhibitor harmine at concentrations of 0 (DMSO), 10 and 20 μM for 48 hours. (b) Images showing the keloid or normal fibroblasts after 48 hours of treatment. Bar = 100 μm (c) Volcano plot showing the differentially expressed genes between 10 and 0 μM treatment obtained by bulk RNA-seq analysis for keloid (left panel) or normal (right panel) fibroblasts. Two replicates were prepared for each treatment. q-value < 0.05, the absolute of log2 fold change > 1. (d) Functional enrichment analysis showing the downregulated pathways between 10 and 0 μM treatment in keloid fibroblasts. (e) Functional enrichment analysis showing the downregulated pathways between 10 and 0 μM treatment in normal fibroblasts. (f) The effects of harmine on the mRNA expression of TWIST1 and genes encoding representative fibrotic factors in keloid fibroblasts. (g) The effects of harmine on the mRNA expression of TWIST1 and genes encoding representative fibrotic factors in normal fibroblasts. In f and g, ∗q < 0.05, ∗∗q < 0.01, q-values were calculated based on differentially expression analysis of RNA-seq data. (h) Immunoblot assays of protein extract from keloid fibroblasts following harmine treatment (left panel) and densitometry analysis of the immunoblot data (right panel). (i) Immunoblot assays of protein extract from normal fibroblasts following harmine treatment (left panel) and densitometry analysis of the immunoblot data (right panel). In h and i, Control means the group without any treatment. Data represent mean ± SEM of three biological replicates. ∗P < 0.05, ∗∗ P < 0.01. One-way ANOVA with LSD post hoc tests. ns: not significant; qRT-PCR, quantitative real-time reverse transcriptase–PCR; RNA-seq, RNA-sequencing.

       Cell-cell communication analysis reveals ligand-receptor interaction changes in keloid fibroblasts

      To define the cell-cell communication landscape and identify its alterations in keloids, we performed analysis using CellPhoneDB 2.0 (
      • Efremova M.
      • Vento-Tormo M.
      • Teichmann S.A.
      • Vento-Tormo R.
      CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit receptor-ligand complexes.
      ). A dense communication network among fibroblasts, VECs, neural cells and keratinocytes was identified, which constituted a core signaling module in the skin. Compared with the normal condition (Supplementary Figure S5a left panel), the total number of interactions for most lineages increased in keloids (Supplementary Figure S5a right panel), reflecting enhanced intercellular communications in diseased conditions. In both conditions, the most abundant interactions occurred between fibroblasts and VECs (Supplementary Figure S5b). Furthermore, we identified the ligand-receptor pairs showing significant changes in specificity between any of the nonfibroblast lineages and fibroblasts in keloid versus normal conditions (Supplementary Table S11). As shown in Supplementary Figure S5c, the significantly altered signaling pathways included a wealth of fibrosis-related signals, such as TGFB1, NOTCH1 and FGF signaling. Notably, TGFB1, a key ligand of the TGFβ signaling pathway implicated in keloid formation (
      • Peltonen J.
      • Hsiao L.L.
      • Jaakkola S.
      • Sollberg S.
      • Aumailley M.
      • Timpl R.
      • et al.
      Activation of collagen gene expression in keloids: co-localization of type I and VI collagen and transforming growth factor-beta1 mRNA.
      ), was found to be secreted by lymphatic endothelial cells, leukocytes, mural cells, neural cells, smooth muscle cells and VECs in keloids. The TGFB1-TGFβ receptor 2 interactions between these cell lineages and fibroblasts became significantly more specific in keloids compared with normal conditions (permutation test P-value < 0.05). Also, we explored the alterations of ligand signals broadcast by fibroblasts (Supplementary Figure S5d). Notably, fibroblasts may affect other cells in keloids through alterations in ligand-receptor interactions involved in NOTCH signaling.

       Cellular heterogeneity and regulatory changes in keloid VECs

      We also examined the heterogeneity and regulatory changes in keloid VECs. Unbiased clustering revealed four clusters with distinct expression profiles (Supplementary Figure S6a and Supplementary Table S12), which could be marked with specific marker combinations: c2-PECAM1+ HMOX1high CXCL3-, c4-PECAM1+ ACKR1high HMOX1-, c5-PECAM1+ CXCL12high CXCL3- and c18-PECAM1+ CXCL3+ (Supplementary Figure S6b). The presence of the four subpopulations was confirmed by immunofluorescence staining (Supplementary Figure S7). Correlation analysis revealed that c4 appeared to be distant from the other clusters (Supplementary Figure S6c). Functional enrichment analysis revealed that c4 represented antigen-presenting endothelial cells expressing major histocompatibility complex class II markers (Supplementary Figure S6d). Then, we identified the lineage-specific differentially regulated pathways of VECs in keloid versus normal skin tissue (Supplementary Figures S6e and S8 and Supplementary Table S13). As expected, the VEGFR signaling pathway, which is important in pathological angiogenesis (
      • Shibuya M.
      Vascular endothelial growth factor (VEGF) and its receptor (VEGFR) signaling in angiogenesis: a crucial target for anti- and pro-angiogenic therapies.
      ), was significantly activated. Intriguingly, tumor-related signaling pathways, such as oncogenic MAPK signaling, signaling by WNT in cancer and the regulation of PTEN gene transcription were activated in VECs in keloids. Furthermore, we identified the ligand-receptor pairs that were altered in keloids between VECs and the other cell lineages (Supplementary Figure S6f and g). Notably, EFNB2-EPHA4, a ligand-receptor pair of Eph-ephrin signaling that is known for promoting sprouting angiogenesis (
      • Kania A.
      • Klein R.
      Mechanisms of ephrin-Eph signalling in development, physiology and disease.
      ), was significantly altered in keloids. Also, keloids VECs may regulate the transcriptional states of fibroblasts and smooth muscle cells through EFNB2-EPHA4 signaling (Supplementary Figure S6g).

      Discussion

      Understanding the cellular heterogeneity and regulatory changes in tissues under diseased conditions is fundamental to successful medical therapy development. Our single-cell analysis provided a wealth of candidate molecules for developing targeted approaches for combating the excessive ECM deposition caused by fibroblasts, or for suppressing the active angiogenesis mediated by dysregulated VECs in keloids. Substantial evidence indicated that the TGFβ signaling pathway plays a central role in keloids or other fibrotic disorders (
      • Piersma B.
      • Bank R.A.
      • Boersema M.
      Signaling in fibrosis: TGF-β, WNT, and YAP/TAZ converge.
      ). The TGFB1 ligand can reprogram fibroblasts to myofibroblasts, which synthesize large volumes of collagen fibers, a hallmark of keloids (
      • Ashcroft K.J.
      • Syed F.
      • Bayat A.
      Site-specific keloid fibroblasts alter the behaviour of normal skin and normal scar fibroblasts through paracrine signalling.
      ;
      • Nangole F.W.
      • Agak G.W.
      Keloid pathophysiology: fibroblast or inflammatory disorders?.
      ). TWIST1 functions as a profibrotic factor in a TGFβ/SMAD3/p38-dependent manner (
      • Ning X.
      • Zhang K.
      • Wu Q.
      • Liu M.
      • Sun S.
      Emerging role of Twist1 in fibrotic diseases.
      ;
      • Palumbo-Zerr K.
      • Soare A.
      • Zerr P.
      • Liebl A.
      • Mancuso R.
      • Tomcik M.
      • et al.
      Composition of TWIST1 dimers regulates fibroblast activation and tissue fibrosis.
      ). In addition, TGFβ signaling has been implicated in tumor angiogenesis (
      • Pardali E.
      • ten Dijke P.
      Transforming growth factor-beta signaling and tumor angiogenesis.
      ). In this study, we observed the activation of TGFβ receptor signaling in both fibroblasts (Figure 3b) and VECs (Supplementary Figure S6e) from keloid tissue. SMAD3, the key downstream regulator of TGFβ signaling, was found to be critical for the diseased states of keloid fibroblasts through independent analyses (Figure 5d and f). The TGFB1-TGFβ receptor 2 interactions between nonfibroblast lineages and fibroblasts were found to become significantly more specific in keloids (Supplementary Figure S5c). Together, these results consolidated our view that TGFβ signaling could serve as a medical target for simultaneously suppressing myofibrogenesis and angiogenesis in keloids.
      Similarly, Eph-ephrin signaling represents another promising target pathway for regulating the states of keloid fibroblasts and VECs. We observed the activation of Eph-ephrin signaling in both fibroblasts (Figure 3b) and VECs (Supplementary Figure S6e) from keloid tissue. We found that EPHB2, a receptor involved in Eph-ephrin signaling, was ranked at the top of the network of keloid fibroblasts (Figure 5e). The angiogenesis-promoting EFNB2-EPHA4 interactions (
      • Kania A.
      • Klein R.
      Mechanisms of ephrin-Eph signalling in development, physiology and disease.
      ) between VECs and other cells were significantly altered in keloids (Supplementary Figure S6f and g).
      Although keloids are generally regarded as benign dermal tumors, keloids display malignant features (
      • Lim C.P.
      • Phan T.T.
      • Lim I.J.
      • Cao X.
      Stat3 contributes to keloid pathogenesis via promoting collagen production, cell proliferation and migration.
      ). The mechanism underlying these features has not been fully elucidated. The PTEN gene encodes a negative regulator of the PI3K/ACT/mTOR pathway that controls cell proliferation and survival (
      • Chalhoub N.
      • Baker S.J.
      PTEN and the PI3-kinase pathway in cancer.
      ). The dysfunction of PTEN (a known tumor suppressor) has been frequently observed in malignant tumors and some fibrotic disorders, such as lung fibrosis (
      • Tian Y.
      • Li H.
      • Qiu T.
      • Dai J.
      • Zhang Y.
      • Chen J.
      • et al.
      Loss of PTEN induces lung fibrosis via alveolar epithelial cell senescence depending on NF-κB activation.
      ). Genetic studies revealed that polymorphisms in PTEN were associated with an increased risk of keloids in the Han Chinese population (
      • Li J.
      • Chen X.
      • Zhou Y.
      • Xie H.
      • Li J.
      • Su J.
      • et al.
      Risk of keloid associated with polymorphic PTEN haplotypes in the Chinese Han population.
      ). Reduced expression of PTEN has been observed in keloid tissue (
      • Sang P.F.
      • Wang H.
      • Wang M.
      • Hu C.
      • Zhang J.S.
      • Li X.J.
      • et al.
      NEDD4-1 and PTEN expression in keloid scarring.
      ). We found that the pathway of the negative regulation of PTEN transcription was activated in both keloid fibroblasts (Figure 3b) and VECs (Supplementary Figure S6e), which indicates that the dysregulation of the PTEN-mediated pathway may be responsible for the malignant features of keloids. In addition, other tumor-related signaling pathways were activated in keloid VECs (Supplementary Figure S6e), which supports the view that overlap exists in the dysregulated pathways between keloid and malignant tumors. Our findings have implications for the clinical treatment of keloids: medical therapies applied in tumors, such as drugs targeting MAPK signaling in cancer (
      • Lee S.
      • Rauch J.
      • Kolch W.
      Targeting MAPK signaling in cancer: mechanisms of drug resistance and sensitivity.
      ), may also be effective in keloid treatment.

      Materials & Methods

      More Materials and Methods can be found online as supplementary information.

       Ethical approval

      All human patient recruitment and tissue sampling procedures complied with the ethical regulations approved by Peking Union Medical College Hospital (#ZS-2162). Written informed consent was received from each subject.

       Patients and sample preparation

      Keloid lesional tissues were harvested during plastic surgery from four patients confirmed to show clinical evidence of keloids. No patient had received chemotherapy or radiotherapy prior to surgery. As a control, matched relatively normal skin tissues adjacent to the keloid scar were also sampled.

       Data availability statement

      Datasets related to this article can be found at https://ngdc.cncb.ac.cn/gsa-human/browse/HRA000425, hosted at Genome Sequence Archive.

      ORCIDs

      Conflict of Interest

      The authors state no conflict of interest.

      Acknowledgments

      This work was supported by grants from the Natural Science Foundation of China (81870229, 81900282).

      Author Contributions

      Conceptualization: XLi; Data Curation: XLi; Formal Analysis: XLi; Funding Acquisition: ZZ; Investigation: XLi, WC, QZ, BM; Methodology: XLi, WC, QZ, BM; Project Administration: ZZ, XLo; Resources: ZL, TM, JC, NY; Software: XLi; Supervision: ZZ; Validation: WC; Visualization: XLi; Writing - Original Draft Preparation: XLi, WC; Writing - Review and Editing: XLi, WC, XLo, ZZ

      Supplementary Material

      Supplementary Materials and Methods

       Tissue dissociation

      Enzymatic digestion and mechanical cutting were performed to dissociate the skin tissue into single-cell suspensions. In brief, the skin tissue was washed twice with DMEM. After removal of the adipose tissue under the reticular dermis, the samples were minced into small pieces. Multiplex enzyme was used to digest the samples overnight in 37 °C water bath according to the manufactory’s protocol (GentleMACS Human Whole Skin dissociation kit, 130-101-540, Miltenyi Biotec, Germany). Next, mechanical cutting was performed on a gentleMACS Octo Dissociator using the h_skin_01 program. Thereafter, fragments and large clumps were removed by filtering through a 100-μm filter (BD Falcon) and then a 40-μm filter. A Dead Cell Removal Kit (130-090-101, Miltenyi Biotec) was used to remove cells with low viability. The live cells were resuspended in PBS containing 0.04% BSA. Subsequently, the live cells were centrifuged at 300  relative centrifugal force (rcf) for 5  minutes at 4  °C to obtain a cell pellet. The cell pellet was resuspended and washed twice before diluted to 1 × 10E6 cells per millimeter. All the above steps were completed under aseptic conditions.

       Single-cell RNA-sequencing library preparation and sequencing

      Single-cell Gel Beads-in-Emulsion generation, barcoding, post- Gel Beads-in-Emulsion-RT cleanup, cDNA amplification and cDNA library construction were performed using Chromium Single Cell 3’ Reagent Kit v3 chemistry (10× Genomics, Pleasanton, CA) following the manufacturer’s protocol. The resulting libraries were sequenced in a NovaSeq 6000 System (Illumina, San Diego, CA).

       Sample demultiplexing, barcode processing and unique molecular identifier counting

      The official Cell Ranger software version 3.0.2 (https://support.10xgenomics.com) was applied for sample demultiplexing, barcode processing and unique molecular identifier (UMI) counting. Briefly, the raw base call files generated by the sequencers were demultiplexed into reads in FASTQ format using the cellranger mkfastq pipeline. Then, the reads were processed using the cellranger count pipeline to generate a gene-barcode matrix for each library. During this step, the reads were aligned to the human reference genome (version: GRCh38). The resulting gene-cell UMI count matrices of all samples were ultimately concatenated into a single matrix using the cellranger aggr pipeline.

       Data cleaning, normalization, feature selection, integration and scaling

      The concatenated gene-cell barcode matrix was imported into Seurat version 3.1.0 (
      • Butler A.
      • Hoffman P.
      • Smibert P.
      • Papalexi E.
      • Satija R.
      Integrating single-cell transcriptomic data across different conditions, technologies, and species.
      ;
      • Stuart T.
      • Butler A.
      • Hoffman P.
      • Hafemeister C.
      • Papalexi E.
      • Mauck 3rd, W.M.
      • et al.
      Comprehensive integration of single-cell data.
      ) for data preprocessing. To exclude genes that were likely detected from random noise, we filtered out genes exhibiting counts in fewer than 3 cells. To exclude poor-quality cells that might have resulted from doublets or other technical noise, we filtered cell outliers (> third quartile + 1.5 × interquartile range or < first quartile – 1.5 × interquartile range) based on the number of expressed genes, the sum of UMI counts and the proportion of mitochondrial genes. To further remove doublets, we filtered out cells based on the predictions provided by Scrublet (
      • Wolock S.L.
      • Lopez R.
      • Klein A.M.
      Scrublet: computational identification of cell doublets in single-Cell transcriptomic data.
      ). In addition, cells that were enriched in hemoglobin gene expression were considered red blood cells and were excluded from further analyses. The sum of the UMI counts for each cell was normalized to 10,000 and log-transformed. For each sample, 2,000 features (genes) were selected using the FindVariableFeatures function of Seurat under the default settings. To correct for potential batch effects and identify shared cell states across datasets, we integrated all the datasets via canonical correlation analysis implemented in Seurat. To mitigate the effects of uninteresting sources of variation (e.g., cell cycle), we regressed out the mitochondrial gene proportion, UMI count, S-phase score and G2M phase score (calculated by the CellCycleScoring function) with linear models using the ScaleData function. Then, the data were centered for each gene by subtracting the average expression of that gene across all cells, and were scaled by dividing the centered expression by the SD.

       Dimensional reduction and clustering

      The expression of the selected genes was subjected to linear dimensional reduction through principal component analysis. Then, the first 30 principal components of the principal component analysis were used to compute a neighborhood graph of the cells. The neighborhood graph was ultimately embedded in two-dimensional space using the nonlinear dimensional reduction method of uniform manifold approximation and projection (
      • Becht E.
      • McInnes L.
      • Healy J.
      • Dutertre C.A.
      • Kwok I.W.H.
      • Ng L.G.
      • et al.
      Dimensionality reduction for visualizing single-cell data using UMAP.
      ). The neighborhood graph of cells was clustered using Louvain clustering (resolution = 0.6) (
      • Blondel V.D.
      • Guillaume J.L.
      • Lambiotte R.
      • Lefebvre E.
      Fast unfolding of communities in large networks.
      ).

       Differential expression and functional enrichment analysis

      Differentially expressed genes between two groups of cells were detected with the likelihood-ratio test (test.use: bimod) implemented in the FindMarkers function of Seurat. The significance threshold was set to an adjusted P-value < 0.05 and a log2-fold change > 0.25. Functional enrichment analyses of a list of genes were performed using ClueGO (
      • Bindea G.
      • Mlecnik B.
      • Hackl H.
      • Charoentong P.
      • Tosolini M.
      • Kirilovsky A.
      • et al.
      ClueGO: a cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks.
      ) with an adjusted P-value threshold of 0.05.

       Gene set enrichment analysis

      All the expressed genes were pre-ranked by Signal2Noise (the difference of means between CASE and CTRL scaled by the SD). Then, the ranked gene list was imported into the GSEA software (version: 4.0.1) (
      • Subramanian A.
      • Tamayo P.
      • Mootha V.K.
      • Mukherjee S.
      • Ebert B.L.
      • Gillette M.A.
      • et al.
      Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.
      ). An false discovery rate q-value < 0.05 was considered to be statistically significant. Precompiled gene sets including REACTOME pathways and GENE ONTOLOGY biological processes in MSigDB (version: 7.0) (
      • Liberzon A.
      • Birger C.
      • Thorvaldsdóttir H.
      • Ghandi M.
      • Mesirov J.P.
      • Tamayo P.
      The molecular signatures database (MSigDB) hallmark gene set collection.
      ) were used in this analysis. The results were visualized using the EnrichmentMap plugin of Cytoscape.

       Pseudo-temporal ordering of single cells

      The pseudo-temporal ordering of the cells along the differentiation trajectory was performed using Monocle2 (
      • Qiu X.
      • Hill A.
      • Packer J.
      • Lin D.
      • Ma Y.A.
      • Trapnell C.
      Single-cell mRNA quantification and differential analysis with Census.
      ). Briefly, the ordering was based on 1,000 genes showing differences in expression between clusters selected via the unsupervised dpFeature procedure. Then, the data space was reduced to two dimensions with the DDRTree method. The cells were ultimately ordered in pseudotime, and cells exhibiting high expression of myofibroblast markers were considered to represent the end of the trajectory.
      To identify significantly branch-dependent genes in their expression, we used branched expression analysis modeling, and the statistically significance threshold was set to a q-value < 1E-04.

       Gene regulatory network analysis based on single-cell transcriptomes

      Gene regulatory networks were constructed from single-cell datasets and compared using the method implemented in bigScale2 (
      • Iacono G.
      • Massoni-Badosa R.
      • Heyn H.
      Single-cell transcriptomics unveils gene regulatory network plasticity.
      ). Briefly, gene regulatory networks for CASE and CTRL were inferred with the compute.network function (clustering = direct, quantile. P = 0.90) separately. Genes encoding ribosomal proteins or mitochondrial proteins were excluded from this analysis. Then, the number of edges was homogenized throughout the obtained networks using the homogenize.networks function. Finally, changes in node centralities (the relative importance of genes in the network) in CASE compared to CTRL were identified using the compare.centrality function. Four measures of centrality: degree, betweenness, closeness and pagerank, were considered. The networks were ultimately visualized with Cytoscape (version: 3.7.0).

       Cell-cell communication analysis based on single-cell transcriptomes

      To analyze cell-cell communication based on single-cell transcriptomic datasets, we used CellPhoneDB 2.0 (
      • Efremova M.
      • Vento-Tormo M.
      • Teichmann S.A.
      • Vento-Tormo R.
      CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes.
      ), which contains a curated repository of ligand-receptor interactions and a statistical framework for inferring lineage-specific interactions. Briefly, potential ligand-receptor interactions were established based on the expression of a receptor by one lineage and a ligand by another. Only ligands and receptors expressed in greater than 10% of the cells in any given lineage were considered. The labels of all cells were randomly permuted 1,000 times and the means of the average ligand-receptor expression in the interacting lineages were calculated, thus generating a null distribution for each ligand-receptor pair in each pairwise comparison between lineages. Ultimately, a P-value for the likelihood of lineage specificity for a given ligand-receptor pair was obtained.

       Isolation, culture and treatment of skin fibroblasts

      The dishes were coated by 0.2% gelation for at least 30 minutes at 37 °C. After epidermal removal, skin tissue was minced into small pieces. Then, the pieces were dissolved in DMEM (11965-084, Gibco, Waltham, MA) containing 20% fetal bovine serum (10500064, Gibco) and were seeded onto gelation-coated culture dishes until attachment at 37 °C. The fibroblasts before passage 7 were used in the subsequent experiments. When the confluence of cells reached >90%, the medium was replaced with culture medium with DMSO or harmine (SH8340, Solarbio, China) at different concentrations (10 and 20 μM). After treatment for 48 hours, the cells were harvested for total RNA extraction using TRIzol reagent (15596026, Invitrogen, Waltham, MA). The total RNA was used for bulk RNA-sequencing library construction using a NEBNext Ultra RNA Library Prep Kit for Illumina (#E7530L, NEB, San Diego, CA). Sequencing was performed on an Illumina X Ten system.

       Bulk RNA-sequencing data analysis

      After quality control using fastp (
      • Chen S.
      • Zhou Y.
      • Chen Y.
      • Gu J.
      fastp: an ultra-fast all-in-one FASTQ preprocessor.
      ), sequencing reads were subjected to transcript abundance quantification with kallisto (
      • Bray N.L.
      • Pimentel H.
      • Melsted P.
      • Pachter L.
      Near-optimal probabilistic RNA-seq quantification [published correction appears in Nat Biotechnol 2016;34:888].
      ). The identification of differentially expressed genes between groups was performed using sleuth (
      • Pimentel H.
      • Bray N.L.
      • Puente S.
      • Melsted P.
      • Pachter L.
      Differential analysis of RNA-seq incorporating quantification uncertainty.
      ). The statistical significance threshold was set to a q value of 0.05, and the biological significance threshold was set to a fold change of ±2-fold.

       Multiplexed immunofluorescence staining

      Multiplexed immunofluorescence staining was performed by following previously described methods (
      • Liu S.L.
      • Bian L.J.
      • Liu Z.X.
      • Chen Q.Y.
      • Sun X.S.
      • Sun R.
      • et al.
      Development and validation of the immune signature to predict distant metastasis in patients with nasopharyngeal carcinoma.
      ). Briefly, PANO 7-plex IHC Kits (cat 0004100100, Panovue, Beijing, China) were used according to the manufacturer’s instructions. Following incubation with primary antibody, skin sections were incubated with horseradish peroxidase-conjugated secondary antibody. Then, the sections were subjected to tyramide signal amplification using TSA Fluorescence Kits (Panovue, Beijing, China). For each additional marker, the protocol was repeated by adding an antigen retrieval step. Finally, all the sections were stained for DAPI (D9542, Sigma-Aldrich, St.Louis, MO). The stained slides were scanned using the Mantra System (PerkinElmer, Waltham, MA). All the antibodies used are as follows: Periostin (POSTN, Abcam, Cambridge, United, Kingdom, ab14041), ABCA8 (Invitrogen, Waltham, MA, PA5-60866), c-Rel (REL, Invitrogen, MA5-15859), RGS2 (Invitrogen, PA5-28162), WISP2 (Invitrogen, PA5-13232), CXCL3 (Invitrogen, PA5-103136), CXCL12 (Invitrogen, MA5-23759), HMOX1 (Invitrogen, MA1-112), DARC (ACKR1, Invitrogen, PA5-20360), CD31 (PECAM1, Abcam, ab28364) and PDGFRA (Cell Signaling Technology, Danvers, MA, 3174S).

       Quantitative real-time reverse transcriptase–PCR

      Keloid fibroblasts were treated with 0, 10 and 20 μM harmine for 48 hours. Total RNA was extracted from fibroblasts using the TRIzol reagent (15596018, Invitrogen). To quantify the relative expression of target genes, quantitative real-time reverse transcriptase–PCR with SYBR green as reporter dye was performed following the methods described previously (
      • Xiang Q.
      • Xu F.
      • Li Y.
      • Liu X.
      • Chen Q.
      • Huang J.
      • et al.
      Transcriptome analysis and functional identification of adipose-derived mesenchymal stem cells in secondary lymphedema.
      ). The housekeeping gene GAPDH was used as a control. Three biological replicates (keloid fibroblasts from three patients) were set for each treatment.

       Immunoblot assay

      Whole cell extracts of fibroblasts were prepared using radioimmunoprecipitation assay buffer (P0013B, Beyotime, Shanghai, China) supplemented with protease and phosphatase inhibitor cocktail (P1046, Beyotime). Protein concentrations were quantified by BCA assay (23227, Thermo Fisher, Waltham, MA). Lysates were diluted in ×5 SDS-PAGE sample loading buffer and boiled at 95 °C for 10 minutes. The proteins were separated by 10% SDS-polyacrylamide gel electrophoresis or 4–12% NuPAGE Bis-Tris SDS-PAGE gel electrophoresis (NP0322BOX, Invitrogen), and then were transferred onto NC membranes. Membranes were blocked in 3% BSA (ST023, Beyotime) or 5% nonfat powdered milk (P0216, Beyotime) in ×1 tris-buffered saline (ST673, Beyotime) for 1 hour at room temperature followed by overnight incubation at 4 °C with primary antibody diluted in primary antibody dilution buffer (P0256, Beyotime) and 1 hour incubation with HRP-conjugated secondary antibody diluted at 1:10000 in secondary antibody dilution buffer (P0258, Beyotime) followed by detection with BeyoECL Star (P0018A, Beyotime). The following primary antibodies were used: GAPDH (1:6000, Cell Signaling Technology Cat# 2118L), Collagen I (1:4000, Abcam Cat# ab138492), Collagen III (1:2000, Abcam Cat# ab184993), Smad2 /Smad3 (1:1000, Cell Signaling Technology Cat# 5339S,), p-Smad2 /Smad3 (1:500, Cell Signaling Technology Cat# 8828S) and TWIST1 (1:100, Abcam Cat# ab50887). Densitometry analysis was performed by quantifying the intensity of the bands using ImageJ software.
      Figure thumbnail fx2
      Supplementary Figure S1UMAP plots showing the distribution of the cells from each sample in each cluster and each fibroblast subpopulation. (a) Distribution of the cells from each sample in each cluster. (b) Distribution of the cells from each fibroblast subpopulation. Each sample was subsampled for equal number of cells. The color codes are in accordance with those used in b or c. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; UMAP, uniform manifold approximation and projection.
      Figure thumbnail fx3
      Supplementary Figure S2Differential proportional analysis of each cell lineage in CASE versus CTRL. A significance threshold of a P-value < 0.05 for paired t-tests was used. The relative proportion of each lineage was calculated in each sample. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FB, fibroblast; lEndo, lymphatic endothelial cell; vEndo, vascular endothelial cell; vSMC, vascular smooth muscle cell.
      Figure thumbnail fx4
      Supplementary Figure S3Immunofluorescence staining of the five subpopulations of keloid fibroblasts in keloid tissue. (a) Immunofluorescence staining of the s0 subpopulation. (b) Immunofluorescence staining of the s1 subpopulation. (c) Immunofluorescence staining of the s2 subpopulation. (d) Immunofluorescence staining of the s3 subpopulation. (e) Immunofluorescence staining of the s4 subpopulation. Low-magnification overviews of representative tissue sections together with high-magnification views of the selected regions are shown. Green boxes indicate the selected regions. Bar (long) = 2 mm; Bar (intermediate) = 20 μm; Bar (short) = 5 μm.
      Figure thumbnail fx5
      Supplementary Figure S4The expression of fibrosis-related genes in keloid fibroblasts following the treatment with harmine. Quantitative real-time reverse transcriptase–PCR confirmed that TWIST1 inhibitor harmine could significantly suppresses the expression of fibrosis-related genes in keloid fibroblasts from three patients. Keloid fibroblasts were treated with TWIST1 inhibitor harmine at concentrations of 0, 10 and 20 μM for 48 hours. Paired t-test was used. ∗: P-value adjusted for multiple tests < 0.05, ∗∗: adjusted P-value < 0.01. ID, identification; P, patient.
      Figure thumbnail fx6
      Supplementary Figure S5Cell-cell communications among the cell lineages in keloids and relatively normal skin tissues. (a) Inter-lineage communication networks in keloids (CASE; right panel) and relatively normal skin tissues (CTRL; left panel). The total number of communications is shown for each cell lineage. The line color indicates that the ligands that are broadcast by the cell lineage in the same color. The line thickness is proportional to the number of broadcast ligands. (b) Heatmap showing the number of communications between any two lineages in CASE (right panel) and CTRL (left panel). (c) The ligand-receptor pairs showing significant changes in specificity between any of the nonfibroblast lineages and fibroblasts in CASE versus CTRL. Fibroblasts express receptors and receive ligand signals from other lineages. The dot size reflects the P-value of the permutation tests for lineage-specificity. The dot color denotes the mean of the average ligand-receptor expression in the interacting lineages. (d) The ligand-receptor pairs shown significant changes in specificity between fibroblasts and any one of the nonfibroblast lineages in CASE versus CTRL. Fibroblasts express ligands and broadcast ligand signals for other lineages. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FB, fibroblast; KTR, keratinocyte; lEndo, lymphatic endothelial cell; LEU, leukocyte; MLA, melanocyte; SGC, sweat gland cell; SMC, smooth muscle cell; vEndo, vascular endothelial cell.
      Figure thumbnail fx7
      Supplementary Figure S6The heterogeneity and regulatory changes in vascular endothelial cells from keloid skin tissues. (a) Heatmap showing distinct expression profiles among the four vascular endothelial cell clusters. (b) Representative molecular signatures of the four cell clusters. (c) Correlations among the four cell clusters. (d) Functional enrichment for each of the four clusters. (e) Enrichment plots (upper panel) and leading-edge gene expression heatmaps (the top 20 genes; lower panel) for representative signaling pathways upregulated in vascular endothelial cells in keloid tissue. NES used to compare the analysis results across gene sets. The vertical lines in the enrichment plot show where the members of the gene set appear in the ranked list of genes. Leading-edge genes: the subset of genes in the gene set that contribute most to the observed enrichment. The average expression across the cells in each group is shown in the heatmap. (f) The ligand-receptor pairs showing significant changes in specificity between any of the nonvascular endothelial lineages and vascular endothelial cells in CASE versus CTRL. Vascular endothelial cells express receptors and receive ligand signals from other lineages. The dot size reflects the P-value of the permutation tests for lineage-specificity. The dot color denotes the mean ligand-receptor expression in the interacting lineages. (g) The ligand-receptor pairs showing significant changes in specificity between vascular endothelial cells and any of the nonvascular endothelial lineages in CASE versus CTRL. Vascular endothelial cells express ligands and broadcast ligand signals for other lineages. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FB, fibroblast; KTR, keratinocyte; lEndo, lymphatic endothelial cell; LEU, leukocyte; MLA, melanocyte; NES, normalized enrichment score; SGC, sweat gland cell; SMC, smooth muscle cell; vEndo, vascular endothelial cell.
      Figure thumbnail fx8
      Supplementary Figure S7Immunofluorescence staining of the four subpopulations of vascular endothelial cells in keloid skin tissues. Arrows indicate the target cells. Green boxes indicate the selected regions. Low-magnification overviews of representative tissue sections together with high-magnification views of the selected regions are shown. Bar (long) = 2 mm; Bar (short) = 20 μm.
      Figure thumbnail fx9
      Supplementary Figure S8Gene set enrichment analysis reveals lineage-specific dysregulated pathways for vascular endothelial cells in keloid versus normal skin tissue. The size of the dot reflects the size of the gene set. The dots in red denote upregulated pathways, and dots in blue represent downregulated pathways. A significance threshold of an FDR q-value of 0.05 was used. CASE, keloid lesional skin tissue; CTRL, adjacent normal tissue; FDR, false discovery rate.

      References

        • Al-Tamari H.M.
        • Dabral S.
        • Schmall A.
        • Sarvari P.
        • Ruppert C.
        • Paik J.
        • et al.
        FoxO3 an important player in fibrogenesis and therapeutic target for idiopathic pulmonary fibrosis.
        EMBO Mol Med. 2018; 10: 276-293
        • Ashcroft K.J.
        • Syed F.
        • Bayat A.
        Site-specific keloid fibroblasts alter the behaviour of normal skin and normal scar fibroblasts through paracrine signalling.
        PLoS One. 2013; 8: e75600
        • Chalhoub N.
        • Baker S.J.
        PTEN and the PI3-kinase pathway in cancer.
        Annu Rev Pathol. 2009; 4: 127-150
        • Darby I.A.
        • Laverdet B.
        • Bonté F.
        • Desmoulière A.
        Fibroblasts and myofibroblasts in wound healing.
        Clin Cosmet Investig Dermatol. 2014; 7: 301-311
        • Der E.
        • Suryawanshi H.
        • Morozov P.
        • Kustagi M.
        • Goilav B.
        • Ranabathou S.
        • et al.
        Tubular cell and keratinocyte single-cell transcriptomics applied to lupus nephritis reveal type I IFN and fibrosis relevant pathways [published correction appears in Nat Immunol 2019;20:1556].
        Nat Immunol. 2019; 20: 915-927
        • Efremova M.
        • Vento-Tormo M.
        • Teichmann S.A.
        • Vento-Tormo R.
        CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes.
        Nat Protoc. 2020; 15: 1484-1506
        • Gauglitz G.G.
        • Korting H.C.
        • Pavicic T.
        • Ruzicka T.
        • Jeschke M.G.
        Hypertrophic scarring and keloids: pathomechanisms and current and emerging treatment strategies.
        Mol Med. 2011; 17: 113-125
        • Glass 2nd, D.A.
        Current understanding of the genetic causes of keloid formation.
        J Investig Dermatol Symp Proc. 2017; 18: S50-S53
        • Guerrero-Juarez C.F.
        • Dedhia P.H.
        • Jin S.
        • Ruiz-Vega R.
        • Ma D.
        • Liu Y.
        • et al.
        Single-cell analysis reveals fibroblast heterogeneity and myeloid-derived adipocyte progenitors in murine skin wounds.
        Nat Commun. 2019; 10: 650
        • He H.
        • Suryawanshi H.
        • Morozov P.
        • Gay-Mimbrera J.
        • Del Duca E.
        • Kim H.J.
        • et al.
        Single-cell transcriptome analysis of human skin identifies novel fibroblast subpopulation and enrichment of immune subsets in atopic dermatitis.
        J Allergy Clin Immunol. 2020; 145: 1615-1628
        • Holm A.
        • Heumann T.
        • Augustin H.G.
        Microvascular mural cell organotypic heterogeneity and functional plasticity.
        Trends Cell Biol. 2018; 28: 302-316
        • Iacono G.
        • Massoni-Badosa R.
        • Heyn H.
        Single-cell transcriptomics unveils gene regulatory network plasticity.
        Genome Biol. 2019; 20: 110
        • Johnson N.C.
        • Dillard M.E.
        • Baluk P.
        • McDonald D.M.
        • Harvey N.L.
        • Frase S.L.
        • et al.
        Lymphatic endothelial cell identity is reversible and its maintenance requires Prox1 activity.
        Genes Dev. 2008; 22: 3282-3291
        • Joost S.
        • Zeisel A.
        • Jacob T.
        • Sun X.
        • La Manno G.
        • Lönnerberg P.
        • et al.
        Single-cell transcriptomics reveals that differentiation and spatial signatures shape epidermal and hair follicle heterogeneity.
        Cell Syst. 2016; 3: 221-237.e9
        • Kalucka J.
        • de Rooij L.P.M.H.
        • Goveia J.
        • Rohlenova K.
        • Dumas S.J.
        • Meta E.
        • et al.
        Single-cell transcriptome atlas of murine endothelial cells.
        Cell. 2020; 180: 764-779.e20
        • Kania A.
        • Klein R.
        Mechanisms of ephrin-Eph signalling in development, physiology and disease.
        Nat Rev Mol Cell Biol. 2016; 17: 240-256
        • Kulkarni A.
        • Anderson A.G.
        • Merullo D.P.
        • Konopka G.
        Beyond bulk: a review of single cell transcriptomics methodologies and applications.
        Curr Opin Biotechnol. 2019; 58: 129-136
        • Lee S.
        • Rauch J.
        • Kolch W.
        Targeting MAPK signaling in cancer: mechanisms of drug resistance and sensitivity.
        Int J Mol Sci. 2020; 21: 1102
        • Li J.
        • Chen X.
        • Zhou Y.
        • Xie H.
        • Li J.
        • Su J.
        • et al.
        Risk of keloid associated with polymorphic PTEN haplotypes in the Chinese Han population.
        Wounds. 2014; 26: 21-27
        • Liang X.
        • Ma L.
        • Long X.
        • Wang X.
        LncRNA expression profiles and validation in keloid and normal skin tissue.
        Int J Oncol. 2015; 47: 1829-1838
        • Lim C.P.
        • Phan T.T.
        • Lim I.J.
        • Cao X.
        Stat3 contributes to keloid pathogenesis via promoting collagen production, cell proliferation and migration.
        Oncogene. 2006; 25: 5416-5425
        • Liu X.
        • Chen W.
        • Li W.
        • Li Y.
        • Priest J.R.
        • Zhou B.
        • et al.
        Single-cell RNA-seq of the developing cardiac outflow tract reveals convergent development of the vascular smooth muscle cells.
        Cell Rep. 2019; 28: 1346-1361.e4
        • Mari W.
        • Alsabri S.G.
        • Tabal N.
        • Younes S.
        • Sherif A.
        • Simman R.
        Novel insights on understanding of keloid scar: article review.
        J Am Coll Clin Wound Spec. 2016; 7: 1-7
        • Miller A.J.
        • Du J.
        • Rowan S.
        • Hershey C.L.
        • Widlund H.R.
        • Fisher D.E.
        Transcriptional regulation of the melanoma prognostic marker melastatin (TRPM1) by MITF in melanocytes and melanoma.
        Cancer Res. 2004; 64: 509-516
        • Nakashima M.
        • Chung S.
        • Takahashi A.
        • Kamatani N.
        • Kawaguchi T.
        • Tsunoda T.
        • et al.
        A genome-wide association study identifies four susceptibility loci for keloid in the Japanese population.
        Nat Genet. 2010; 42: 768-771
        • Nangole F.W.
        • Agak G.W.
        Keloid pathophysiology: fibroblast or inflammatory disorders?.
        JPRAS Open. 2019; 22: 44-54
        • Ning X.
        • Zhang K.
        • Wu Q.
        • Liu M.
        • Sun S.
        Emerging role of Twist1 in fibrotic diseases.
        J Cell Mol Med. 2018; 22: 1383-1391
        • Onoufriadis A.
        • Hsu C.K.
        • Ainali C.
        • Ung C.Y.
        • Rashidghamat E.
        • Yang H.S.
        • et al.
        Time series integrative analysis of RNA sequencing and microRNA expression data reveals key biologic wound healing pathways in keloid-prone individuals.
        J Invest Dermatol. 2018; 138: 2690-2693
        • Palumbo-Zerr K.
        • Soare A.
        • Zerr P.
        • Liebl A.
        • Mancuso R.
        • Tomcik M.
        • et al.
        Composition of TWIST1 dimers regulates fibroblast activation and tissue fibrosis.
        Ann Rheum Dis. 2017; 76: 244-251
        • Pardali E.
        • ten Dijke P.
        Transforming growth factor-beta signaling and tumor angiogenesis.
        Front Biosci (Landmark Ed). 2009; 14: 4848-4861
        • Peltonen J.
        • Hsiao L.L.
        • Jaakkola S.
        • Sollberg S.
        • Aumailley M.
        • Timpl R.
        • et al.
        Activation of collagen gene expression in keloids: co-localization of type I and VI collagen and transforming growth factor-beta1 mRNA.
        J Invest Dermatol. 1991; 97: 240-248
        • Piersma B.
        • Bank R.A.
        • Boersema M.
        Signaling in fibrosis: TGF-β, WNT, and YAP/TAZ converge.
        Front Med (Lausanne). 2015; 2: 59
        • Qiu X.
        • Hill A.
        • Packer J.
        • Lin D.
        • Ma Y.A.
        • Trapnell C.
        Single-cell mRNA quantification and differential analysis with Census.
        Nat Methods. 2017; 14: 309-315
        • Reyfman P.A.
        • Walter J.M.
        • Joshi N.
        • Anekalla K.R.
        • McQuattie-Pimentel A.C.
        • Chiu S.
        • et al.
        Single-cell transcriptomic analysis of human lung provides insights into the pathobiology of pulmonary fibrosis.
        Am J Respir Crit Care Med. 2019; 199: 1517-1536
        • Sang P.F.
        • Wang H.
        • Wang M.
        • Hu C.
        • Zhang J.S.
        • Li X.J.
        • et al.
        NEDD4-1 and PTEN expression in keloid scarring.
        Genet Mol Res. 2015; 14: 13467-13475
        • Shibuya M.
        Vascular endothelial growth factor (VEGF) and its receptor (VEGFR) signaling in angiogenesis: a crucial target for anti- and pro-angiogenic therapies.
        Genes Cancer. 2011; 2: 1097-1105
        • Shih B.
        • Garside E.
        • McGrouther D.A.
        • Bayat A.
        Molecular dissection of abnormal wound healing processes resulting in keloid disease.
        Wound Repair Regen. 2010; 18: 139-153
        • Su S.A.
        • Yang D.
        • Wu Y.
        • Xie Y.
        • Zhu W.
        • Cai Z.
        • et al.
        EphrinB2 regulates cardiac fibrosis through modulating the interaction of Stat3 and TGF-β/Smad3 signaling.
        Circ Res. 2017; 121: 617-627
        • Tian Y.
        • Li H.
        • Qiu T.
        • Dai J.
        • Zhang Y.
        • Chen J.
        • et al.
        Loss of PTEN induces lung fibrosis via alveolar epithelial cell senescence depending on NF-κB activation.
        Aging Cell. 2019; 18: e12858
        • Wang J.
        • Wu H.
        • Xiao Z.
        • Dong X.
        Expression profiles of lncRNAs and circRNAs in keloid.
        Plast Reconstr Surg Glob Open. 2019; 7: e2265
        • Yochum Z.A.
        • Cades J.
        • Mazzacurati L.
        • Neumann N.M.
        • Khetarpal S.K.
        • Chatterjee S.
        • et al.
        A first-in-class twist1 inhibitor with activity in oncogene-driven lung cancer.
        Mol Cancer Res. 2017; 15: 1764-1776
        • Zheng G.X.
        • Terry J.M.
        • Belgrader P.
        • Ryvkin P.
        • Bent Z.W.
        • Wilson R.
        • et al.
        Massively parallel digital transcriptional profiling of single cells.
        Nat Commun. 2017; 8: 14049