Shunuo Zhang
Yixin Zhang*
Peiru Min- Department of Plastic and Reconstructive Surgery, Shanghai Ninth People’s Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, Shanghai, China
Hypertrophic scar (HS) is a common skin disorder characterized by excessive extracellular matrix (ECM) deposition. However, it is still unclear how the cellular composition, cell-cell communications, and crucial transcriptionally regulatory network were changed in HS. In the present study, we found that FB-1, which was identified a major type of fibroblast and had the characteristics of myofibroblast, was significantly expanded in HS by integrative analysis of the single-cell and bulk RNA sequencing (RNA-seq) data. Moreover, the proportion of KC-2, which might be a differentiated type of keratinocyte (KC), was reduced in HS. To decipher the intercellular signaling, we conducted the cell-cell communication analysis between the cell types, and found the autocrine signaling of HB-1 through COL1A1/2-CD44 and CD99-CD99 and the intercellular contacts between FB-1/FB-5 and KC-2 through COL1A1/COL1A2/COL6A1/COL6A2-SDC4. Almost all the ligands and receptors involved in the autocrine signaling of HB-1 were upregulated in HS by both scRNA-seq and bulk RNA-seq data. In contrast, the receptor of KC-2, SDC4, which could bind to multiple ligands, was downregulated in HS, suggesting that the reduced proportion of KC-2 and apoptotic phenotype of KC-2 might be associated with the downregulation of SDC4. Furthermore, we also investigated the transcriptionally regulatory network involved in HS formation. The integrative analysis of the scRNA-seq and bulk RNA-seq data identified CREB3L1 and TWIST2 as the critical TFs involved in the myofibroblast of HS. In summary, the integrative analysis of the single-cell RNA sequencing (scRNA-seq) and bulk RNA-seq data greatly improved our understanding of the biological characteristics during the HS formation.
Introduction
Hypertrophic scar (HS) caused by pathologically excessive collagen deposition from skin fibroblasts in the dermis and subcutis are major types of pathological scars that can be regarded as complications to abnormal wound healing (Haverstock, 2001). Unlike keloids, HS does not extend beyond the borders of the original wound, as myofibroblasts, the main effector cell contributing to dermal fibrosis, are present in hypertrophic scarring, and α-smooth muscle actin (α-SMA) is expressed in a nodular formation in the fibroblasts from HS (Lee et al., 2004; Feng et al., 2020). Meanwhile, abnormal keratinocyte differentiation and proliferation, and significantly increased acanthosis can be observed along with hypertrophic scarring (Niessen et al., 2001).
Through animal models, it is found that inflammatory cells, bone marrow-derived fibrocytes, and several peptide-related compounds may serve as essential players in hypertrophic scar development (Satoyoshi, 1966; Song et al., 2020; Wang et al., 2020). It has been well-recognized that the expressions of keratinocyte-derived interleukin-1 (IL-1), tumor necrosis factor-a (TNF-a), platelet-derived growth factor (PDGF), transforming growth factor-b (TGF-b) and basic fibroblast growth factor (bFGF) are associated with extracellular matrix (ECM) remodeling, and medicines or potential reagents to prevent and treat hypertrophic scars could result in inflammatory inhibition via targeting these molecules (Limandjaja et al., 2018; Wang et al., 2020). A previous research has described that fibrocytes, functioning as antigen-presenting cells, may contribute to the upregulation of the inflammatory response, and another study has demonstrated that the TGFβ1/Smad pathway could regulate ECM synthesis via stimulating fibroblasts and inducing fibroblast differentiation into myofibroblasts, thus promoting the formation of hypertrophic scars (Yang et al., 2005; Liu et al., 2012; Shirakami et al., 2020). Also, activation of PI3K/AKT pathway could lead to hypertrophic scar formation and ECM deposition via upregulated expression of Collagen I, Collagen III, α-SMA, and Cleaved caspase-3 in hypertrophic scar fibroblasts (Xiao, 2020; Zhi et al., 2021). In addition, some transcription factors (TFs) such as CRBE3L1 and TWIST2 have been identified as critical regulators in skin diseases (Crespo et al., 2021; Deng et al., 2021). Importantly, the two TFs were also involved in kidney fibrosis (Grunz-Borgmann et al., 2017; Yamamoto et al., 2021). However, few hypertrophic scar animal models could perfectly reflect human skin injuries, and though much effort has been dedicated to the study of abnormal wound healing, the pathogenesis of hypertrophic scarring has not been fully unveiled, and breakthrough development in the therapeutic management for hypertrophic scars is urgently needed, as hypertrophic scarring would lead to considerable morbidity (Wang et al., 2011; Domergue et al., 2015).
Utilization of single-cell RNA sequencing (scRNA-seq) has helped examine cell variability across tumors and interactions between detected cell types (Kumar et al., 2018), and with the aid of scRNA-seq technologies, exploring fibroblast heterogeneity at a single-cell resolution and cell-cell communication in the context of hypertrophic scarring has now become a reality. Recent study identifies serine proteases as regulators of myofibroblast differentiation using single cell sequencing technology (Vorstandlechner et al., 2020). In the present study, we hope to identify potential cell-cell communications and abnormal transcriptionally regulatory network in hypertrophic scarring.
Materials and Methods
Data Collection
The scRNA-seq data of normal and HS samples were collected from Gene Expression Omnibus (GEO) with accession number GSE156326 (Vorstandlechner et al., 2020). The bulk RNA sequencing data was generated by this study. The five normal and five HS skin tissue samples were collected from Nineth People’s Hospital of Shanghai Jiao Tong University, School of Medicine, which was approved by the Human Research Ethics Committee of this hospital. The written informed consents were collected from each patient. All samples were stored in −80°C for the following experiments.
RNA Extraction, Library Construction and RNA Sequencing Analysis
RNA extraction and sequencing was performed at the Beijing Genomics Institute (BGI), using an Illumina HiSeq 4,000 sequencer (Illumina) following the standard manufacturer’s protocols as described by the previous study (Hillen et al., 2019). The raw data was first preprocessed by excluding the reads with low quality. The clean reads were the mapped to human reference genome (GRCh37/hg19) using Hisat v2.2.1 (Kim et al., 2015). The gene expression quantification was performed by using Stringtie v2.1.4 (Pertea et al., 2015) and R ballgown package (Frazee et al., 2015). The default parameters of Hisat and Stringtie were used in this study. The count table was generated by the python script prepDE.py (https://ccb.jhu.edu/software/stringtie/dl/prepDE.py) from Stringtie. The R/Bioconductor DESeq2 (Love et al., 2014) package was applied to differential gene expression analysis.
Cell Clustering Analysis
The unique molecular identifiers (UMIs) count-based scRNA-seq data of six human skin samples from GEO accession GSE156326 were used for the cell clustering analysis, which was implemented in R Seurat package. Cells with less than 200 UMIs were eliminated and features detected in less than three cells were filtered. The multiple datasets were integrated by SCTransform in Seurat package. The top 3,000 highly variable features were selected by FindVariableFeatures. The clusters were found at a resolution of 0.8 by FindClusters, and T-distributed Stochastic Neighbor Embedding (t-SNE) was applied to reduce the dimensionality. The cell-type marker genes were detected by FindAllMarkers function at adjusted p-value < 0.05, minimal percentage difference >0.25, and log2 fold change >0.5. All the marker genes of the cell clusters were collected from the earlier study (Vorstandlechner et al., 2020). This analysis was implemented by R Seurat package (Stuart et al., 2019).
Estimation of Cell Proportion
The cell proportions were estimated using MuSiC(Wang et al., 2019), which used a deconvolution method based on marker genes of cell types and gene expression matrices of both scRNA-seq and bulk RNA-seq to estimate the cell proportions of bulk RNA-seq data. The count-based expression data of both scRNA-seq and bulk RNA-seq was applied to this analysis.
Gene Set Enrichment Analysis (GSEA)
The GSEA was conducted against KEGG, Reactome (Jassal et al., 2020), WikiPathway, and Gene Ontology (GO-bp). The GSEA was implemented in R clusterProfiler (Yu et al., 2012).
Cell-Cell Communication Analysis
The cell-cell communication was predicted by R CellChat (Jin et al., 2021) package. The normalized count and cell types by Seurat were used for this analysis. This analysis was separately conducted on the normal and HS cells.
Transcriptionally Regulatory Network Analysis
The transcriptionally regulatory network analysis was performed by R SCENIC package (Aibar et al., 2017). The AUCell values for the transcription factors (TFs) were used for differential analysis (Aibar et al., 2017).
Statistical Analyses
All the statistical analyses were performed in R (version 4.1.0). The two-sample and pairwise comparisons were conducted by Wilcoxon rank sum test. The symbols of *, **, and *** represent the statistical significance at 0.05, 0.01, and 0.001, respectively.
Results
Single-Cell and Bulk RNA-Seq Reveal the Cellular Diversity and Heterogeneity of Skin Tissues.
To reveal the cellular diversity and heterogeneity of skin tissues, we collected a single-cell RNA sequencing (RNA-seq) dataset of three normal (NS) and three hypertrophic scar (HS) samples from Gene Expression Omnibus (GEO) database (See Materials and methods). After excluding the cells with low quality, we identified 18 cell clusters from 15,276 cells by the dimensional reduction and clustering analysis (Figure 1A), and characterized the cell clusters using the marker genes from previous studies (Solé-Boldo et al., 2020; Vorstandlechner et al., 2020; Liu et al., 2021). The cell clusters included six clusters of fibroblasts (FB-1/2/3/4/5/6, represented by marker genes: FBLN1, COL1A1, APOE, and APCDD1), three clusters of keratinocytes (KC-1/2/3: KRT1 and KRT14), two clusters of endothelial cells (EC-1/2: THBD and SELE), T cells (TC: PTPRC), smooth muscle cells (SMC: ACTA2 and RGS5), Langerhans cells (LA: AIF1), lymphatic endothelial cells (LEC: PROX1 and PDPN), dendritic cells (DC: HLA−DRB1 and CD1C), macrophage (MP: CD68), and melanocyte (MC: MITF) (Figure 1B), suggesting that the cell types were well characterized by the marker genes.
FIGURE 1. The cell populations and marker genes in the skin samples of hypertrophic scar (HS) and normal skins. (A) The cell clusters visualized by the dimensional reduction of t-distributed stochastic neighbor embedding (t-SNE). Each point represents one cell, and the filled colors represent the cell types. (B) The expression specificity of cell type-specific marker genes within the skin samples using scRNA-seq data. (C) The differential expression of the cell type-specific marker genes between HS and normal skins using bulk RNA-seq data.
To reveal the expression patterns of the cell-type marker genes in skin tissues, we also collected five NS and five HS tissue samples for bulk RNA sequencing. As shown in Figure 1C, the representative marker genes were found to be differentially expressed between the NS and HS tissues (adjusted p-value < 0.05 and fold change >1.5). Moreover, about 26.83% of the cell-type marker genes (376/1,401) were observed to be differentially expressed between NS and HS samples. The deregulation of the cell-type marker genes in HS samples indicated that the cell proportions and gene expression patterns were changed in HS.
Differential Proportion Analysis Reveals Significant Expansion of Fibroblast and Contraction of KC Subpopulations in HS
We next attempted to identify HS-associated cell lineages or clusters that were significantly expanded or contracted in HS. The FB-1 accounted for the largest number of cells, followed by FB-2, TC, and EC-2 based on the single-cell RNA-seq data (Figure 2A). Furthermore, we also conducted gene set enrichment analysis on the marker genes of those cell types to characterize their functionalities (Supplementary Table S1). Remarkably, the three major cell types, including FB-2, EC-2, and T cell, were characterized by collagen degradation, VEGFA-VEGFR2 Signaling Pathway, and TYROBP Causal Network, respectively. The cell proportions of the three NS and three HS samples could be estimated by the single-cell RNA-seq data. To improve the reliability of the differential proportion analysis, we also estimated the cell proportions for the ten bulk RNA-seq samples using MuSiC (Wang et al., 2019), a deconvolution method to estimate cell type proportions from bulk RNA sequencing data. The differential proportion analysis revealed that FB-1 and FB-2 were expanded in HS, while the proportion of KC-2 was reduced in HS (Figure 2B). Accordingly, the FB-1 and KC-2 marker genes were enriched in the up- and down-regulated genes in HS by bulk RNA-seq data (Figure 2C), respectively. However, FB-2 marker genes were enriched in neither the upregulated nor the downregulated genes in HS. These results suggested that the changed proportions of these cell types might be closely associated with HS.
FIGURE 2. The cell proportions within HS and normal skin samples. (A) The cell proportions within the integrated HS and normal skin samples. (B) The cell types with differential proportions between HS and normal skin. (C) The differential expression levels of the FB-1 and KC-2 specific marker genes between the HS and normal skin samples.
Fibroblast Heterogeneity in Skin Tissues
As fibroblast clusters were identified as the major cell types in skin tissues by scRNA-seq data, we then investigated their functional differences by comparing their expression profiles. The differential gene expression analysis revealed that FB-1, FB-3, FB-4, and FB-5 had their top-ten specific marker genes, while no marker genes specifically upregulated in FB-2 or FB-6 were identified (Figure 3A). Consistent with the excessive ECM deposition observed in bulk RNA-seq data of HS (Figure 3B), ECM-related pathways, such as extracellular matrix organization, elastic fiber formation, collagen formation, and matrix metalloproteinases, and PI3K-Akt signaling were significantly upregulated in FB-1 (Figure 3C). FB-3 was characterized by the inflammatory pathways, such as photodynamic therapy-induced NF-kB survival signaling, NOD-like receptor signaling pathway, and signaling by interleukins (Figure 3C). The oxidative stress and Wnt signaling were identified as the signature pathway of FB-4 and FB-5 (Figure 3C), respectively. These results indicated that fibroblasts were functionally heterogenous in skin tissues.
FIGURE 3. The fibroblast subpopulations in the skin samples. (A) The expression specificity of the FB subpopulation specific genes across the FB cell populations. (B) The upregulated pathways in FB-1 cells of HS. (C) The subpopulation specific biological functions across the FB subpopulations.
As FB-1 was one of the major cell types expanded in HS, we then investigated the biological function of this FB type. Interestingly, the FB-1 cells from HS samples could be clearly distinguished from those from NS samples (Figure 4A). The graph-based clustering analysis revealed that the C1 subpopulation was significantly enriched by the cells from HS (Figure 4A, hypergeometric test, p-value < 0.05). Moreover, we also found that the myofibroblast gene signatures were higher in C1 and HS, as compared with C2 and NS cells (Figure 4B), respectively, suggesting that the FB-1 from HS might have a myofibroblast phenotype. Moreover, the myofibroblast signature genes, including COL12A1, COL1A1, COL3A1, CTHRC1, and PCSK1N, were expressed higher in C1 than C2 (Figure 4C). The bulk RNA-seq data further confirmed the upregulation of the five genes except PCSK1N (Figure 4D). Notably, the five genes were involved in collagen formation, which is a critical regulator of myofibroblast differentiation (Yuen et al., 2010). These results indicated that HS formation was associated with both the increased proportion and the biological characteristics of myofibroblast.
FIGURE 4. Biological function of HS-related FB-1 cluster. (A) The clustering of the HS and normal skin (NS)-related FB-1 cells grouped by disease (HS vs NS) and graph-based classification (C1 vs C2), respectively. (B) The differential module score of myofibroblast cells between HS and NS, and C1 and C2. (C) The expression levels of myofibroblast cell marker genes including COL12A1, COL1A1, COL3A1, CTHRC1, and PCSK1N in the HB-1 cells. (D) The differential expression levels of the myofibroblast cell marker genes between HS and NS bulk samples.
Dysregulation of Keratinocytes Signature Genes in HS
As the proportion of type 2 keratinocyte was reduced in HS, we then asked whether KC-2 in HS showed distinct expression patterns relative to those in NS. Compared to KC-1 and KC-3, the KC-2 exerted distinct expression profiles (Figure 5A). The gene set enrichment analysis revealed that the three KC cell types (KC-1/2/3) might be associated with antigen presentation, keratinization, and transforming growth factor-beta (TGFbeta) signaling (Figure 5B). The tSNE analysis revealed that the KC-2 in HS and NS could not be clearly distinguished (Figure 5C), suggesting that the KC-2 in HS and NS showed similar expression patterns. The differential gene expression analysis between KC-2 cells of HS and NS revealed that cell senescence or apoptosis-related pathways such as aging, neuron death, and regulation of neuron apoptotic process were upregulated in KC-2 cells of HS (Figure 5D), suggesting that the keratinocyte apoptosis might be associated with the reduced proportion of KC-2 in HS.
FIGURE 5. The keratinocyte subpopulations in the skin samples. (A) The expression specificity of the KC subpopulation specific genes across the KC cell populations. (B) The subpopulation specific biological functions across the KC subpopulations. (C) The clustering of the KC-2 cells grouped by disease (HS vs NS). (D) The representative biological functions in KC-2 cells of HS.
Cell-Cell Communications
To decipher the intercellular signaling, we conducted the cell-cell communication analysis between the cell types using CellChat (Jin et al., 2021). As the FB1-1 and KC-2 exhibited dominant functional roles in HS formation, we then investigated the signaling sources of the two cell types. Specifically, we identified 13 ligand-receptor (LR) pairs upregulated in FB-1 of HS and 17 LR pairs downregulated in KC-2 of HS (Figure 6A, adjusted p-value < 0.05). Particularly, we found that the autocrine signaling of HB-1 through COL1A1/2-CD44 and CD99-CD99 achieved the high communication probabilities (Figure 6A). Similarly, the cell-cell contacts between FB-1/FB-5 and KC-2 through COL1A1/COL1A2/COL6A1/COL6A2-SDC4 had the high communication probabilities in the downregulated cell-cell communication network (Figure 6A). The expression analysis revealed that the ligands and receptors such as COL1A1, COL1A2, CD44, and CD99 involved in the signaling transduction of FB-1 had higher expression levels in FB-1, as compared with the other cell types (Figure 6B). In contrast, the receptors involved in the signaling transduction of KC-2 were specifically expressed in KC-2 and KC-3, and the ligands were more specifically expressed in FB-1 and FB-5 (Figure 6B). Moreover, we also investigated the differential expression levels of those ligands and receptors between the HS and NS using the bulk RNA-seq data. Notably, all the ligands and receptors involved in signaling transduction of FB-1 except SELE were found to be upregulated in HS (Figure 6C). Interestingly, the receptor SDC4, which could bind to multiple ligands, were downregulated in HS, suggesting that the reduced proportion of KC-2 and apoptotic phenotype of KC-2 might be associated with the downregulation of SDC4.
FIGURE 6. The differential cell-cell communications between HS and NS. (A) The upregulated and downregulated signaling in HS. The x- and y-axis represent the cell-cell communications and ligand-receptor pairs, respectively. (B) The expression levels of the ligands and receptors involved in the upregulated or downregulated signaling across the cell types. (C) The differential expression levels of the ligands and receptors between the HS and NS bulk RNA-seq samples.
Transcriptionally Regulatory Network Involved in HS Formation
To infer the transcriptionally regulatory network involved in HS formation, we estimated the transcriptional activities of the transcription factors (TFs) using SCENIC(Aibar et al., 2017). The differential transcriptional activity analysis revealed that the FB-1 cells in HS had higher transcriptional activities of CREB3L1 and TWIST2 than those in NS (Figure 7A). Unfortunately, we did not observe any TFs which showed different activities between HS and NS of KC-2. HS had a higher proportion of the FB-1 cells expressing CREB3L1 than NS as CREB3L1 was lowly expressed in FB-1 (Figure 7B). The FB-1 in HS had a higher expression level of TWIST2 than that in NS (Figure 7B, adjusted p-value < 0.05). The previous study identified JUN as a critical regulator by promoting hypertrophic skin scarring (Griffin et al., 2021). Consistently, we also found that the transcriptional activity of JUN was higher in HS than NS (Supplementary Figure S1).
FIGURE 7. The transcription factors (TFs) and regulatory network involved in HS. (A) The differential transcriptional activities of the key TFs in FB-1 of HS samples. (B) The differential proportion of FB-1 cells expressing CREB3L1 and differential expression of TWIST2 between the FB1- cells of HS and NS. (C) The differential expression of TFs between the bulk RNA-seq samples of HS and NS. (D,E) Upregulation of TF target genes in the bulk RNA-seq samples of HS by gene set enrichment analysis (GSEA).
Furthermore, we also investigated the expression levels of the TFs and target genes in HS and NS samples using bulk RNA-seq data. Consistently, the two TFs, CREB3L1 and TWIST2, were upregulated in HS as compared with NS using the bulk RNA-seq data (Figure 7C, p-value < 0.01). In according with the finding using scRNA-seq data, the target genes of CREB3L1 or TWIST2 also gathered in the upregulated genes in HS (Figures 7D,E, FDR <0.05). Similarly, the target genes of the two TFs were also significantly overlapped with the upregulated genes in the HS samples from a previous study (GEO accession: GSE188952, Supplementary Figure S2), further indicating that the transcriptional activities of the two TFs were upregulated in HS. These results suggested that the two TFs, CREB3L1 and TWIST2, might play key roles in HS formation.
Discussion
Hypertrophic scar (HS) is a fibroproliferative skin disorder characterized by excessive extracellular matrix (ECM) deposition. However, it is still unclear how the cellular composition, cell-cell communications, and crucial transcriptionally regulatory network were changed in HS. In the present study, we found that fibroblasts (FB) and keratinocytes (KC) had multiple cell clusters. The FB clusters were characterized by excessive ECM deposition, inflammatory, and proliferative signatures. Particularly, excessive ECM deposition, identified as the signature of FB-1, is associated with myofibroblast (Hinz and Lagares, 2020) and fibrosis (Agarwal et al., 2013). Moreover, the FB-1 was expanded in HS and FB-1 accounted for the largest proportion in the HS samples, suggesting that an increase of myofibroblast proportion in HS might be the major cause of the HS. The myofibroblast was considered as a therapeutic target of HS (Feng et al., 2020; Yang et al., 2020). The three KC cell types (KC-1/2/3) might be associated with antigen presentation, keratinization, and transforming growth factor-beta (TGF-beta) signaling. Remarkably, the proportion of KC-2, which might be a differentiated type of KC, was reduced in HS. The cell senescence or apoptosis-related pathways such as aging, neuron death, and regulation of neuron apoptotic process were upregulated in KC-2 cells of HS, suggesting that the keratinocyte apoptosis might be associated with the reduced proportion of KC-2 in HS. It has been well recognized that keratinocyte proliferation and migration are critical for re-epithelialization during cutaneous wound healing (DiPersio et al., 2016; Piipponen et al., 2020).
To decipher the intercellular signaling, we conducted the cell-cell communication analysis between the cell types, and found the autocrine signaling of HB-1 through COL1A1/2-CD44 and CD99-CD99 and the intercellular contacts between FB-1/FB-5 and KC-2 through COL1A1/COL1A2/COL6A1/COL6A2-SDC4. It has been reported that COL1A1/2 expression promotes the myofibroblast differentiation (Gambini et al., 2018), fibrosis (Kong et al., 2019), and HS formation (Bi et al., 2017). The receptor, CD44, was associated with an abnormal accumulation of extracellular matrix (Messadi and Bertolami, 1993). Moreover, CD44 was also identified as a mediator of myofibroblast differentiation and fibrosis (Woods et al., 2021). CD99, the other key receptor expressed in myofibroblast, was implicated in skin disorders (Kazakov et al., 2005; Belonogov and Ruksha, 2013). However, the intracellular signaling regulated by CD99 was still unclear. In contrast, the receptors involved in the signaling transduction of KC-2 were specifically expressed in KC-2 and KC-3, and the ligands were more specifically expressed in FB-1 and FB-5. Interestingly, the receptor SDC4, which could bind to multiple ligands, was downregulated in HS, suggesting that the reduced proportion of KC-2 and apoptotic phenotype of KC-2 might be associated with the downregulation of SDC4. Particularly, SDC4 could induce the keratinocyte activation (Bizzarro et al., 2019; Pessolano et al., 2019), which played a vital role in wound healing (Zhang et al., 2021). The low expression of SDC4 indicated that the signaling from fibroblasts could not be transduced to keratinocyte.
Furthermore, we also investigated the transcriptionally regulatory network involved in HS formation. The integrative analysis of the scRNA-seq and bulk RNA-seq data identified CREB3L1 and TWIST2 as the critical TFs involved in the myofibroblast of HS. Consistently, CRBE3L1 has been identified as one of the critical TFs involved in fibrotic skin diseases by a scRNA-seq study (Deng et al., 2021). The other TF, TWIST2, was relevant to Setleis syndrome (SS), a focal facial dermal dysplasia presenting with bilateral temporal skin lesions (Crespo et al., 2021). Moreover, TWIST1, an important paralog of TWIST2, was involved in the fibrogenesis of keloid fibroblasts, and might serve as a therapeutic target of keloid (Liu et al., 2021), suggesting that TWIST2 might also be considered as a promising therapeutic target in HS.
Even though cell-cell communications and key TFs involved in HS have been predicted based on the integrative analysis, there is still a lack of enough evidence to link the relationship between the transcriptionally regulatory network of CRBE3L1 or TWIST2 and the intercellular signaling like COL1A1/2-CD44 and CD99-CD99, which will be clarified by our future study.
In summary, the present study identified that an increase of myofibroblast proportion in HS, abnormal ligand-receptor interaction, and two key TFs might be promising therapeutic targets for HS. The integrative analysis of the scRNA-seq and bulk RNA-seq data not only identified key cell types altered in HS, but also shed light on the potential cell-cell communications and intracellular TF regulatory network involved in HS, which greatly improved our understanding of the biological characteristics during the HS formation.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.biosino.org/node/project/detail/OEP002674, OEP002674
Ethics Statement
The studies involving human participants were reviewed and approved by the Human Research Ethics Committee of Nineth People’s Hospital of Shanghai Jiao Tong University, School of Medicine. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
PM and YZ designed the study. SZ and YZ collected the samples and conducted data analysis. YZ and PM interpreted the data. SZ wrote the manuscript. PM and YZ revised the manuscript.
Funding
This work was supported by National Natural Science Foundation of China (81801918).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.806740/full#supplementary-material
References
Agarwal, P., Schulz, J.-N., Blumbach, K., Andreasson, K., Heinegård, D., Paulsson, M., et al. (2013). Enhanced Deposition of Cartilage Oligomeric Matrix Protein Is a Common Feature in Fibrotic Skin Pathologies. Matrix Biol. 32 (6), 325–331. doi:10.1016/j.matbio.2013.02.010
Aibar, S., González-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., et al. (2017). SCENIC: Single-Cell Regulatory Network Inference and Clustering. Nat. Methods 14 (11), 1083–1086. doi:10.1038/nmeth.4463
Belonogov, R. N., and Ruksha, T. G. (2013). Characteristics of CD99 Expression in the Skin of Psoriasis Patients. Arkh Patol 75 (5), 26–29.
Bi, S., Chai, L., Yuan, X., Cao, C., and Li, S. (2017). MicroRNA-98 Inhibits the Cell Proliferation of Human Hypertrophic Scar Fibroblasts via Targeting Col1A1. Biol. Res. 50 (1), 22. doi:10.1186/s40659-017-0127-6
Bizzarro, V., Belvedere, R., Pessolano, E., Parente, L., Petrella, F., Perretti, M., et al. (2019). Mesoglycan Induces Keratinocyte Activation by Triggering Syndecan‐4 Pathway and the Formation of the Annexin A1/S100A11 Complex. J. Cel Physiol 234 (11), 20174–20192. doi:10.1002/jcp.28618
Crespo, N. E., Torres-Bracero, A., Renta, J. Y., Desnick, R. J., and Cadilla, C. L. (2021). Expression Profiling Identifies TWIST2 Target Genes in Setleis Syndrome Patient Fibroblast and Lymphoblast Cells. Int. J. Environ. Res. Public Health 18 (4). doi:10.3390/ijerph18041997
Deng, C.-C., Hu, Y.-F., Zhu, D.-H., Cheng, Q., Gu, J.-J., Feng, Q.-L., et al. (2021). Single-cell RNA-Seq Reveals Fibroblast Heterogeneity and Increased Mesenchymal Fibroblasts in Human Fibrotic Skin Diseases. Nat. Commun. 12 (1), 3709. doi:10.1038/s41467-021-24110-y
DiPersio, C. M., Zheng, R., Kenney, J., and Van De Water, L. (2016). Integrin-mediated Regulation of Epidermal Wound Functions. Cell Tissue Res 365 (3), 467–482. doi:10.1007/s00441-016-2446-2
Domergue, S., Jorgensen, C., and Noël, D. (2015). Advances in Research in Animal Models of Burn-Related Hypertrophic Scarring. J. Burn Care Res. 36 (5), e259–e266. doi:10.1097/bcr.0000000000000167
Feng, Y., Wu, J.-J., Sun, Z.-L., Liu, S.-Y., Zou, M.-L., Yuan, Z.-D., et al. (2020). Targeted Apoptosis of Myofibroblasts by Elesclomol Inhibits Hypertrophic Scar Formation. EBioMedicine 54, 102715. doi:10.1016/j.ebiom.2020.102715
Frazee, A. C., Pertea, G., Jaffe, A. E., Langmead, B., Salzberg, S. L., and Leek, J. T. (2015). Ballgown Bridges the gap between Transcriptome Assembly and Expression Analysis. Nat. Biotechnol. 33 (3), 243–246. doi:10.1038/nbt.3172
Gambini, E., Perrucci, G. L., Bassetti, B., Spaltro, G., Campostrini, G., Lionetti, M. C., et al. (2018). Preferential Myofibroblast Differentiation of Cardiac Mesenchymal Progenitor Cells in the Presence of Atrial Fibrillation. Translational Res. 192, 54–67. doi:10.1016/j.trsl.2017.11.003
Griffin, M. F., Borrelli, M. R., Garcia, J. T., Januszyk, M., King, M., Lerbs, T., et al. (2021). JUN Promotes Hypertrophic Skin Scarring via CD36 in Preclinical In Vitro and In Vivo Models. Sci. Transl Med. 13 (609), eabb3312. doi:10.1126/scitranslmed.abb3312
Grunz-Borgmann, E. A., Nichols, L. A., Wang, X., and Parrish, A. R. (2017). Twist2 Is Upregulated in Early Stages of Repair Following Acute Kidney Injury. Int. J. Mol. Sci. 18 (2). doi:10.3390/ijms18020368
Hillen, M. R., Pandit, A., Blokland, S. L. M., Hartgring, S. A. Y., Bekker, C. P. J., van der Heijden, E. H. M., et al. (2019). Plasmacytoid DCs from Patients with Sjögren's Syndrome Are Transcriptionally Primed for Enhanced Pro-inflammatory Cytokine Production. Front. Immunol. 10, 2096. doi:10.3389/fimmu.2019.02096
Hinz, B., and Lagares, D. (2020). Evasion of Apoptosis by Myofibroblasts: a Hallmark of Fibrotic Diseases. Nat. Rev. Rheumatol. 16 (1), 11–31. doi:10.1038/s41584-019-0324-5
Jassal, B., Matthews, L., Viteri, G., Gong, C., Lorente, P., Fabregat, A., et al. (2020). The Reactome Pathway Knowledgebase. Nucleic Acids Res. 48 (D1), D498–D503. doi:10.1093/nar/gkz1031
Jin, S., Guerrero-Juarez, C. F., Zhang, L., Chang, I., Ramos, R., Kuan, C.-H., et al. (2021). Inference and Analysis of Cell-Cell Communication Using CellChat. Nat. Commun. 12 (1), 1088. doi:10.1038/s41467-021-21246-9
Kazakov, D. V., Morrisson, C., Plaza, J. A., Michal, M., and Suster, S. (2005). Sarcoma Arising in Hyaline-Vascular Castleman Disease of Skin and Subcutis. The Am. J. dermatopathology 27 (4), 327–332. doi:10.1097/01.dad.0000171606.55810.86
Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a Fast Spliced Aligner with Low Memory Requirements. Nat. Methods 12 (4), 357–360. doi:10.1038/nmeth.3317
Kong, M., Hong, W., Shao, Y., Lv, F., Fan, Z., Li, P., et al. (2019). Ablation of Serum Response Factor in Hepatic Stellate Cells Attenuates Liver Fibrosis. J. Mol. Med. 97 (11), 1521–1533. doi:10.1007/s00109-019-01831-8
Kumar, M. P., Du, J., Lagoudas, G., Jiao, Y., Sawyer, A., Drummond, D. C., et al. (2018). Analysis of Single-Cell RNA-Seq Identifies Cell-Cell Communication Associated with Tumor Characteristics. Cel Rep. 25 (6), 1458–1468. doi:10.1016/j.celrep.2018.10.047
Lee, J. Y.-Y., Yang, C.-C., Chao, S.-C., and Wong, T.-W. (2004). Histopathological Differential Diagnosis of Keloid and Hypertrophic Scar. Am. J. dermatopathology 26 (5), 379–384. doi:10.1097/00000372-200410000-00006
Limandjaja, G. C., van den Broek, L. J., Breetveld, M., Waaijman, T., Monstrey, S., de Boer, E. M., et al. (2018). Characterization of In Vitro Reconstructed Human Normotrophic, Hypertrophic, and Keloid Scar Models. Tissue Eng. C: Methods 24 (4), 242–253. doi:10.1089/ten.tec.2017.0464
Liu, J., Wang, Y., Pan, Q., Su, Y., Zhang, Z., Han, J., et al. (2012). Wnt/β-catenin Pathway Forms a Negative Feedback Loop during TGF-Β1 Induced Human normal Skin Fibroblast-To-Myofibroblast Transition. J. Dermatol. Sci. 65 (1), 38–49. doi:10.1016/j.jdermsci.2011.09.012
Liu, X., Chen, W., Zeng, Q., Ma, B., Li, Z., Meng, T., et al. (2021). Single-cell RNA-Seq Reveals Lineage-specific Regulatory Changes of Fibroblasts and Vascular Endothelial Cells in Keloids. J. Invest. Dermatol.
Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8
Messadi, D. V., and Bertolami, C. N. (1993). CD44 and Hyaluronan Expression in Human Cutaneous Scar Fibroblasts. Am. J. Pathol. 142 (4), 1041–1049.
Niessen, F. B., Andriessen, M. P., Schalkwijk, J., Visser, L., and Timens, W. (2001). Keratinocyte-derived Growth Factors Play a Role in the Formation of Hypertrophic Scars. J. Pathol. 194 (2), 207–216. doi:10.1002/path.853
Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T.-C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie Enables Improved Reconstruction of a Transcriptome from RNA-Seq Reads. Nat. Biotechnol. 33 (3), 290–295. doi:10.1038/nbt.3122
Pessolano, E., Belvedere, R., Bizzarro, V., Franco, P., Marco, I., Petrella, F., et al. (2019). Annexin A1 Contained in Extracellular Vesicles Promotes the Activation of Keratinocytes by Mesoglycan Effects: An Autocrine Loop through FPRs. Cells 8 (7). doi:10.3390/cells8070753
Piipponen, M., Li, D., and Landén, N. X. (2020). The Immune Functions of Keratinocytes in Skin Wound Healing. Int. J. Mol. Sci. 21 (22). doi:10.3390/ijms21228790
Satoyoshi, E.Honda (1966). Demyelinating Cerebrospinal Meningitis. Nihon Ishikai Zasshi 56 (11), 1183–1194.
Shirakami, E., Yamakawa, S., and Hayashida, K. (2020). Strategies to Prevent Hypertrophic Scar Formation: a Review of Therapeutic Interventions Based on Molecular Evidence. Burns Trauma 8, tkz003. doi:10.1093/burnst/tkz003
Solé-Boldo, L., Raddatz, G., Schütz, S., Mallm, J.-P., Rippe, K., Lonsdorf, A. S., et al. (2020). Single-cell Transcriptomes of the Human Skin Reveal Age-Related Loss of Fibroblast Priming. Commun. Biol. 3 (1), 188. doi:10.1038/s42003-020-0922-4
Song, J., Li, X., and Li, J. (2020). Emerging Evidence for the Roles of Peptide in Hypertrophic Scar. Life Sci. 241, 117174. doi:10.1016/j.lfs.2019.117174
Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M., et al. (2019). Comprehensive Integration of Single-Cell Data. Cell 177 (7), 1888–1902. doi:10.1016/j.cell.2019.05.031
Vorstandlechner, V., Laggner, M., Kalinina, P., Haslik, W., Radtke, C., Shaw, L., et al. (2020). Deciphering the Functional Heterogeneity of Skin Fibroblasts Using Single‐cell RNA Sequencing. FASEB j. 34 (3), 3677–3692. doi:10.1096/fj.201902001rr
Wang, J., Ding, J., Jiao, H., Honardoust, D., Momtazi, M., Shankowsky, H. A., et al. (2011). Human Hypertrophic Scar-like Nude Mouse Model: Characterization of the Molecular and Cellular Biology of the Scar Process. Wound Repair Regen. : official Publ. Wound Healing Soc. [and] Eur. Tissue Repair Soc. 19 (2), 274–285. doi:10.1111/j.1524-475x.2011.00672.x
Wang, X., Park, J., Susztak, K., Zhang, N. R., and Li, M. (2019). Bulk Tissue Cell Type Deconvolution with Multi-Subject Single-Cell Expression Reference. Nat. Commun. 10 (1), 380. doi:10.1038/s41467-018-08023-x
Wang, Z.-C., Zhao, W.-Y., Cao, Y., Liu, Y.-Q., Sun, Q., Shi, P., et al. (2020). The Roles of Inflammation in Keloid and Hypertrophic Scars. Front. Immunol. 11, 603187. doi:10.3389/fimmu.2020.603187
Woods, E. L., Grigorieva, I. V., Midgley, A. C., Brown, C. V. M., Lu, Y. A., Phillips, A. O., et al. (2021). CD147 Mediates the CD44s-dependent Differentiation of Myofibroblasts Driven by Transforming Growth Factor-Β1. J. Biol. Chem. 297, 100987. doi:10.1016/j.jbc.2021.100987
Xiao, Y. (2020). MiR-486-5p Inhibits the Hyperproliferation and Production of Collagen in Hypertrophic Scar Fibroblasts via IGF1/PI3K/AKT Pathway. J. Dermatolog Treat., 1–10. doi:10.1080/09546634.2020.1728210
Yamamoto, A., Morioki, H., Nakae, T., Miyake, Y., Harada, T., Noda, S., et al. (2021). Transcription Factor Old Astrocyte Specifically Induced Substance Is a Novel Regulator of Kidney Fibrosis. FASEB J. : official Publ. Fed. Am. Societies Exp. Biol. 35 (2), e21158. doi:10.1096/fj.202001820r
Yang, L., Scott, P. G., Dodd, C., Medina, A., Jiao, H., Shankowsky, H. A., et al. (2005). Identification of Fibrocytes in Postburn Hypertrophic Scar. Wound Repair Regen. : official Publ. Wound Healing Soc. [and] Eur. Tissue Repair Soc. 13 (4), 398–404. doi:10.1111/j.1067-1927.2005.130407.x
Yang, X., Xiao, Y., Zhong, C., Shu, F., Xiao, S., Zheng, Y., et al. (2020). ABT-263 Reduces Hypertrophic Scars by Targeting Apoptosis of Myofibroblasts. Front. Pharmacol. 11, 615505. doi:10.3389/fphar.2020.615505
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118
Yuen, A., Laschinger, C., Talior, I., Lee, W., Chan, M., Birek, J., et al. (2010). Methylglyoxal-modified Collagen Promotes Myofibroblast Differentiation. Matrix Biol. 29 (6), 537–548. doi:10.1016/j.matbio.2010.04.004
Zhang, B., Wu, Y., Mori, M., and Yoshimura, K. (2021). Adipose-Derived Stem Cell Conditioned Medium and Wound Healing: A Systematic Review. Reviews: Tissue engineering Part B.
Keywords: hypertrophic scar, extracellular matrix, single-cell RNA sequencing, cell-cell communication, myofibroblast, keratinocyte
Citation: Zhang S, Zhang Y and Min P (2022) Single-Cell and Bulk Transcriptome Data Integration Reveals Dysfunctional Cell Types and Aberrantly Expressed Genes in Hypertrophic Scar. Front. Genet. 12:806740. doi: 10.3389/fgene.2021.806740
Received: 01 November 2021; Accepted: 09 December 2021;
Published: 03 January 2022.
Edited by:
Tao Huang, Shanghai Institute of Nutrition and Health (CAS), ChinaReviewed by:
Chengfei Zhang, Nanjing Medical University, ChinaXuechao Wan, Northwestern University, United States
Juliang Qin, East China Normal University, China
Copyright © 2022 Zhang, Zhang and Min. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Peiru Min, YXJ1X3JlbkBtc24uY29t; Yixin Zhang, emhhbmd5aXhpbjY2ODhAMTYzLmNvbQ==