Samples from keloid patients were collected after obtaining informed consent and Institutional Review Board (IRB) approval (IRB number: SMC 2020-03-032). Patients confirmed to have clinical evidence of keloids were analyzed in this study.

Keloid samples were classified as AK and IK based on a comprehensive evaluation that included patient feedback such as itching, pain intensity, and an assessment of persistent growth for at least 6 months before and after tissue collection. AK was identified in patients who showed persistent growth at the biopsy or excision site when observed for more than 6 months, increased in size compared to the original lesion, and reported moderate or greater pain and pruritus intensity (6 or greater on a scale of 0-10). IK, on the other hand, includes patients who report no change in size at the biopsy or excision site and mild to no pain (3 or less on a 0-10 scale) when observed for at least 6 months. This classification system was implemented to capture the dynamic nature of keloid development and to distinguish between active and inactive phases. By incorporating these clinical criteria, we aimed to provide a more detailed understanding of the RNA expression patterns associated with different stages of keloid development. We collected six AK samples, four IK samples, and three MS samples from patients with keloids. Keloid samples were obtained on multiple regions via 4 mm punch biopsy or partial excision to minimize entire lesion damage. In this study, we analyzed keloid samples from chest and back regions only, excluding other areas. These samples were used to investigate the cellular microenvironment of keloid.

Single-cell suspensions of keloid and scar samples were subjected to scRNA-seq using Chromium Single Cell Controller and Single Cell 3′ Library & Gel Bead Kit (10X Genomics, Pleasanton, CA, USA) in accordance with the manufacturer’s protocol.

In detail, single-cell suspensions were generated immediately after resection as previously described (31, 32). Briefly, the keloid tissue was washed in phosphate-buffered saline (PBS). Excised keloids were immersed in physiological saline and immediately transferred to the lab. After the sample was minced into small pieces with a scalpel in a Petri dish, prepared dissociation solution (2 mg/ml Liberase TL in PBS, Sigma Aldrich) was used to digest at 37°C for 60 min. The debris was filtered out with a 70-μm cell strainer (#352340, Corning). The remaining tissue was transferred into a dissociation solution containing 0.25% Trypsin for the second round of dissociation. Subsequently, the live cells were resuspended in PBS-0.04% bovine serum albumin. Cells were pelleted by centrifuging at 500 x g for 5 min. Live cells were counted and confirmed using the LUNA-FL dual fluorescence cell counter (Logos Biosystems). To generate single-cell gel bead emulsions, suspensions of skin cells were immediately injected into the Chromium Controller (10x Genomics) targeting 9,000 cells, respectively. Following the manufacturer’s protocols, the Chromium Single Cell 3’ Reagent V3.1 Chemistry Kit (10x Genomics) was used to construct single-cell libraries. Then, prepared libraries were sequenced on the NovaSeq 6000 Sequencing System (Illumina, USA) targeting 50,000 raw reads per cell.

CellRanger pipeline (10x Genomics) was used to map the sequencing reads to a human reference genome (GRCh38). The raw gene expression matrix generated from the CellRanger pipeline was processed for downstream analyses using the Seurat package version 4.2.0 in R version 4.1.1 software (R Foundation for Statistical Computing, Vienna, Austria). Downstream analysis was performed after cell quality filtering. Cells were included if they expressed greater than 500 unique molecular identifier (UMI) counts, fewer than 6,500 genes, greater than 200 genes, and fewer than 15% of mitochondrial gene expression in UMI counts. Doublet identification was performed using the tool DoubletFinder (v2.0.3) by creating artificial doublets and measuring the proportion of artificial k nearest neighbors (pANN) for each cell using PC distance. The cells were then ranked according to the expected number of doublets.

The 10x sequencing data were mapped to human reference genome (GRCh38; official Cell Ranger reference) for generating the raw gene expression matrix. The Seurat package version 4.2 was used for preprocessing and normalizing gene expression data. Genes that were expressed in fewer than 3 cells were excluded in the expression matrix. After quality control filtering, cells from multiple donors were merged using the standard integration protocol described in Seurat v4. Following the Seurat tutorial, we selected the 3000 most variably expressed genes to identify major cell types in each sample. We used the single-cell data integration method “Harmony” to correct for batch effects in scRNA-seq data sets used in this study before performing downstream analysis.

Cells were clustered using dimensionality reduction and Uniform Manifold Approximation and Projection (UMAP) with Seurat R package. The DEGs for each cluster were calculated to ascribe an initial annotation using the Seurat function “FindAllMarkers” (Wilcoxon’s rank sum test) with default parameters. To build second level clusters, cells belonging to subpopulations were reanalyzed separately; all steps were performed using NormalizeData, FindVariableFeature, ScaleData, RunPCA, FindNeighbours, FindClusters, and RunUMAP, the methods implemented in the Seurat package for each of the major cell types. Then, second-level cluster-based doublet exclusion was performed. After doublet removal, we repeated the above mentioned steps to identify subclusters.

For the gene ontology (GO) analysis, we ranked the differentially expressed genes (DEGs) for the SC-0 and SC-2 clusters. Candidate DEGs were further filtered at P-value <0.001 and average log (fold change) > 2. We used the Metascape web tool (www.metascape.org) to conduct GO analysis. The annotation dB R package org.Hs.eg.db was used to map gene identifier. Gene Set Enrichment Analysis (GSEA) was performed using R package clusterProfiler and gene sets from the Molecular Signature Database (MSigDB, gsea-msigdb.org).

To identify shared marker genes for keloids, differential expression was performed using the “FindMarkers” function implemented in Seurat with default parameters. The marker gene set was filtered to include genes with an adjusted p-value <0.05 and an average log fold change higher than 1.6. We further excluded genes with the ribosome, with mitochondria, and with detectable expression in >40% of normal control clusters from these shared keloid marker genes.

The analysis was conducted using the interaction tool NicheNet to explore cell-cell communication. The definition of expressed genes as described in Puram et al. was used to determine background expressed genes in subpopulations. Briefly, NicheNet infers intercellular communication based on the ligand activities and expression of previously defined downstream targets regulated by these prioritized ligands. The function “Predict_ligand_activities” implemented in NicheNet was applied to infer ligand-target potential scores. We also applied Cellchat to infer possible cell-cell ligand-receptor interactions present within keloid cell type. Using the same clusters included for NicheNet analysis, we confirmed which ligands influence transcription during keloid development. A significant interaction was determined when the mean receptor and ligand expression of the subpopulations were significantly higher than those of subpopulations of cells determined from random permutations.

Tissues were fixed with 4% paraformaldehyde and paraffin sections were prepared. Sections were stained with hematoxylin and eosin (H&E). Toluidine blue was used to stain all mature MCs by binding to serglycin proteoglycans in secretory granules. MCs were assessed by quantifying their densities within images of toluidine blue-stained slides in ten random locations of a tissue section (magnification 20×). MC density was expressed as cells/mm².

All specimens were fixed in formalin and embedded in paraffin blocks; the blocks were then cut into 4-µm sections. All tissue sections were immunohistochemical (IHC) stained for the MC marker tryptase and the nerve fiber marker protein gene product 9.5 (PGP9.5). Tissue sections were incubated simultaneously with two primary antibodies: mouse anti-human MC tryptase monoclonal antibody (working dilution 1:100, Santa Cruz Biotechnology) and mouse anti-human PGP9.5 monoclonal antibody (working dilution 1:100, Bio-Rad). Sections were incubated with anti-mouse multimer labeled with horseradish peroxidase (blue staining) for Tryptase and with anti-mouse multimer labeled with alkaline phosphatase substrate (pink staining) for PGP9.5.

All statistical analyses were performed using R4.1.1 software (R Foundation for Statistical Computing, Vienna, Austria). The Mann-Whitney U test was employed for comparing differences in means using the R program. Chi-square tests were used to compare differences in categorical variables among the three groups. Cell counts from AK, IK, and MS samples obtained via scRNAseq were compared with one another with one-way Analysis Of Variance (ANOVA) with post hoc Tukey’s honest significant difference tests. Statistical significance was accepted at p<0.05 for all comparisons.

Histopathological features of AK, IK, and MS tissues were reviewed and evaluated by dermatologists and a pathologist. Remarkably, H&E tissue sections of AK revealed notable increases in cellular components, which were distributed amidst scattered collagens bundles ( Figure 1 ). In contrast, histological analysis of IK samples showed a distinct presentation characterized by densely packed and thick collagen bundles accompanied by sparse cellularity.

To gain a more comprehensive understanding and confirm histopathological findings within the AK compared with IK, we performed a scRNA-seq of fresh keloids derived from three AK patients and two IK patients. Each keloid was divided into center and peripheral areas for the analysis. MS obtained from three keloid patients was also collected during surgery (n=3) ( Supplementary Tables S1 ). We collected two MS samples each from two patients who had already been sampled for IKs.

After quality control and removal of doublets, we profiled a total of 87,051 single cells (71,617 cells from keloids and 11,161 cells from MS). All profiled cells were then clustered and visualized with Unsupervised Uniform Manifold Approximation and Projection (UMAP) ( Figure 2A ; Supplementary Figure S1 ). We then classified cells into 11 major cell types and annotated cell clusters using classic markers and particular transcriptional signatures: FB (COL3A1), vascular endothelial cells (EC; PECAM and VWF), lymphatic endothelial cells (LEC; LYVE1), myofibroblasts (ACTA2 and TAGLN), lymphocytes (LC; CD3D), mast cells (MC; CTSG ⁺), macrophage (MAC; AIF), keratinocytes (KC; KRT1), proliferating keratinocytes (proliferating KC; COL17A1 and MKI67), melanocytes (MEL; PMEL), and neural cells (NC; S100B) ( Figures 2B, C ; Supplementary File S1 ) (33, 34).

Recent studies have described FB heterogeneity within keloids and suggested this molecular heterogeneity is important for understanding underlying mechanisms of keloid pathophysiology (2, 21, 22). To address changes occurring during the fibrotic process in keloids, we conducted an unsupervised clustering analysis (n=23,151) of an FB subpopulation. Results of the scRNA-seq analysis revealed that FBs were divided into four subpopulations using specific markers from previous studies: secretory-papillary, secretory-reticular, mesenchymal, and pro-inflammatory FB ( Supplementary File S2 ) (2, 22, 35). Next, we compared differences between keloids and matured scars obtained from keloid patients ( Figures 2D, E ). Consistent with previous findings, the proportion of mesenchymal subpopulation was increased in keloids when compared to MS. Upregulated genes identified in mesenchymal FBs were primarily associated with ECM organization and skeletal system development, suggesting their potential roles in keloid pathogenesis.

The scRNA-seq analysis showed that immune profiles differed significantly between keloids and MSs, especially in the MC population. In line with previous studies, MCs specifically expressed GATA2, MS4A2, and CPA3 ( Figure 3A ; Supplementary File S3 ) (22, 35). Interestingly, the proportion of MCs was significantly enriched in keloid samples compared to MS samples. Building upon these findings, our subsequent work focused on investigating the MC population within keloids. Toluidine blue immunohistochemical stain, a known marker for MCs, further supported the enrichment pattern of MCs in keloids. This enriched pattern was particularly dominant in active lesions of keloids ( Figure 3B ). These results were also depicted in the box plot, which illustrated MC proportions as determined by scRNA-seq analyses. Our comprehensive analyses revealed a notable abundance of both LC and MAC within AK samples. Collective enrichment of LC, MAC, and MC suggested inflammatory response and immune system activation that could be observed in AK tissues. On the other hand, FB was significantly enriched in keloids, showing a higher abundance in IK samples (p-value = 0.042, Cochran Armitage trend test; Figure 3C ).

We examined the expression pattern of marker genes associated with activated and degranulated MCs, specifically FCER1G, BTK, and GATA2. We found that these marker genes were also significantly elevated in the MCs population of the keloid group compared to the MC population of the MS group, suggesting a potential association between keloid pathogenesis and stressed state of MCs ( Figure 3D ) (36–38). In addition to an increase in the number of MCs, there was an upregulation of tryptase expression and a downregulation of chymase expression within keloid tissues ( Figure 3E ). The upregulation of serine protease tryptase is particularly intriguing as it has been implicated in various functions, including induction of substance P release, activation of neurokinin 1 receptor, and amplification of inflammation, leading to thermal and mechanical hyperalgesia (39). Conversely, downregulation of chymase in skin might contribute to prolonged survival of inflammatory cytokines and neuropeptides, potentially promoting inflammation and itching in keloids.

Our further analyses revealed that MCs derived from AK and IK showed an upregulation of sphingosine 1-phosphate receptors, specifically S1PR1 and S1PR2, along with an upregulation of NTRK1, which encodes TrkA. These receptors are involved in regulation of pain-related nerve growth factor signaling (39). Comparatively, MCs derived from MSs did not exhibit the same level of upregulation, whereas the highest level of expression was observed in MCs from AK ( Figure 3F ). Collectively, our findings suggest that MCs are implicated in both pain and pruritus through expression of neuropeptide receptors (40–42).

Next, we performed subgroup analysis focusing on the population of NCs expressing S100B. Consistent with previous studies, keloids showed a higher proportion of NCs than MSs ( Figure 4A ; Supplementary File S4 ) (43). Further clustering analysis of the NC population revealed the presence of four distinct subclusters. Interestingly, NC-0 and NC-2 were specifically observed in keloid samples ( Figure 4B ). Differentially Expressed Genes (44) analysis confirmed that neural progenitor cell markers such as NES were highly expressed in these clusters ( Figure 4C ) (43). Gene Ontology (GO) enrichment analyses revealed that clusters NC-0 and NC-2 exhibited upregulation of DEGs that were enriched in specific GO terms, including NABA CORE MATRISOME, extracellular matrix organization, and nervous system development ( Supplementary Figure S2 ). These results could be supported by previous findings showing that NCs participate in signal transduction and act on changes in ECM structure to promote keloid formation (43, 45).

Expression of nmSC (non-myelinating Schwann cell) genes such as L1CAM and EMP1 was observed in all lesions, whereas expression of mSC (myelinating Schwann cell) genes including MBP and MPZ was specific to MS. Interestingly, the expression of genes related to neuronal precursor cell genes, which could directly contribute to tissue regeneration and skin repair, was upregulated in the IK group ( Figure 4D ) (46, 47). These results suggest that NCs might have undergone damage during the formation of keloids and that neural progenitor cells are in the process of proliferating to repair the associated nerve damage concurrently. A more detailed analysis of NCs revealed that AK specifically expressed pain sensation-related genes such as TRPA1 (nociception-related sodium channels), NPPB (itch-related neuropeptide), NPR1 (Nppb receptor), HRH1 (histamine receptor), and HRH4 ( Figure 4E ) (27, 48–51). These data indicate that AK has higher expression levels of presumed unmyelinated, peptidergic nociceptors known to have the potential to induce pruritic responses in humans than IK.

To eliminate the possibility of bias, we analyzed only keloid samples from the chest and back, excluding other sites, and also observed the expression of markers between chest and back keloids, which revealed no site-specific differences. Gender-specific differences were not observed either in the analysis ( Supplementary Figure S3 ).

We investigated interactions of MCs and NCs with various cell types through the R package “cellchat” (52). Comparing keloids and MS, our data showed that cell-cell interactions between FBs and other cell types were more enriched in keloids than in MS ( Supplementary Figure S2 ), consistent with previous studies (45). The number of interactions among MCs, NCs, and FBs was measured. Although the difference was not statistically significant, there was a notable increase in the number of interactions between MCs and NCs in keloids compared to MS ( Supplementary Figure S2 ). These results suggest that MCs and NCs might play a role in the maintenance of keloid activation.

In our analyses of IK lesions, FBs exhibited the highest level of combined outgoing and incoming signals among all cell types examined. However, in AK lesions, enhanced interactions were observed across all cell types, highlighting the heterogeneous features associated with keloid activation ( Figure 5A ). Through analysis of cellchat, we observed significant signaling changes that could potentially drive the onset of disease pathogenesis. Specifically, in AK, ligand VEGFA and its receptor VEGFR2 exhibited high activities, particularly in the signaling from MCs to ECs. On the contrary, in IK, we discovered active interactions of TGFB1-(TGFBR1+TGFBR2) from MCs to FBs ( Figure 5B ). We observed the expression of TPSAB1, which encodes Tryptase, in MCs of both AK and IK. Expression of F2RL1, which encodes PAR2, showed significant upregulation in the mesenchymal subpopulation of FBs, particularly in IK ( Figure 5C ). To explore ligand-receptor interactions in more detail, we used ligand-receptor interaction tool NicheNet (53). Using NicheNet, we explored interactions between NC and MC that were significantly enriched in AK ( Figures 5D, E ). The predicted ligands exhibiting the highly expressed in NC from AK included IL33, CCL28, and CXCL14. Through ligand activity analysis, we identified receptors known to be targeted by these ligands in MC, such as NTRK1, IL1RL1, and S1PR1. Furthermore, we mapped genes that were differentially expressed in MC of AK, revealing predicted ligands with elevated expression in MC of AK, including TNFSF14, IL13, CTSG, CXCL2, PTGS2, and LIF. In parallel, we mapped the expression patterns of receptors known to be targeted by these ligands in NC and confirmed the expression of IL13RA1 and IL4R in NC. To support the potential role of these cell interactions, we used PGP 9.5 as a marker for NCs and tryptase for MCs. Through IHC staining with tryptase and PGP 9.5 antibodies, we confirmed this cell-cell interaction, demonstrating the close apposition of MCs and nerve terminals in AKs ( Figure 5F ).