Characteristics of Tumor Immune Gene and Immune Cell Inltration during Keloid Formation

Background: Keloids are benign broproliferative skin tumors that can cause disgurement and disability. Although current research has sought to examine keloids from the perspectives of genetics, inammation, immunity, and tumorigenesis, their pathological mechanisms remain unclear. Methods: In this study, we used three datasets of tumor immune gene expression proling from the normal skin tissue of keloid patients (N group), inammation tissue of keloid patients (I group), and keloid tissue of keloid patients (K group) to describe the occurrence and characteristics of keloid development. Tumor immune-related genes were analyzed, and the differentially expressed genes (DEGs) between the three groups were compared. Gene Ontology (GO) categories and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were carried out to determine the main functions of the differentially expressed genes and keloid-related pathways. Results: We identied several genes that may play an important role in keloid development. These genes are CCR1, SELL, CCR7, CD40LG, CD69, CXCL8, ITGAM, ITGAX, CD86, and CXCL9. GO analysis revealed that there were variations in biological processes (BP) between I group and N group, including regulation of lymphocyte activation and T-cell activation. Similar variations were also found between I group and K group, which may play an important role in keloid initiation and formation. Variations in molecular function (MF) were markedly enriched in cytokine receptor binding and receptor ligand activity. Analysis of the KEGG pathway between I group and N group revealed that DEGs were primarily enriched in cytokine−cytokine receptor interaction and viral protein interaction with cytokine and cytokine receptor. We identied a higher proportion of M2 macrophages in N group than in I group, although the difference was not obvious. M1 macrophage production differed signicantly between I group and K group. The proportions of CD8 + T cells varied signicantly between N group and K group. We traced multiple tumor immune-related hub genes from keloid formation and analyzed immune cell subsets in keloid development. Possible molecular mechanisms were described in this study using bioinformatics.

methods, including surgery, radiotherapy, steroid injection, and other approaches. However, some of these methods, such as surgery and radiotherapy, have caused large trauma to patients. All of the treatments have certain keloid recurrence rates, which make them far from ideal treatment methods. Therefore, it is necessary to further elucidate the pathogenesis of keloid development in order to lay a foundation for clinical exploration of new therapeutic methods.
Elucidating the pathogenesis of keloid development is indispensable to the study of pathogenic genes.
The genetic theory of keloids is mainly founded on the discovery of keloid families. Bayat reported that increased familial clustering in keloids, their increased prevalence in certain races, and increased concordance in identical twins suggest a strong heritable predisposition to keloid formation [4]. Heritable factors for keloid formation may be related to single nucleotide polymorphisms. A genome-wide association study showed that four SNP sites were signi cantly involved in keloid formation. rs8032158SNP is located in intron 5 of the NEDD4 gene on chromosome 15, which leads to abnormal cell proliferation [5]. However, a majority of keloid patients do not have any family history [6,7]. The characteristics of the pathogenic genes of these patients have not been clearly elucidated. The change of pathogenic keloid genes requires further exploration.
It has been reported that the density of CD14 + macrophages and CD3 + cells in keloid tissue is much higher than that of normal skin tissue [8]. The distribution of CD1α + Langerhans cells, CD68 + macrophages, and CD20 + B lymphocytes in keloid tissues are also higher than those of normal skin [9]. Moreover, macrophages in keloid tissues are under stress, and can promote Treg cells by regulating Foxp3 differentiation [8]. These suggest that keloid occurrence and development are closely related to the role of immune cell interaction. However, exploration of tumor immune-related genes and cells in keloid pathogenesis remains poorly understood.
In this study, based on the tumor immune gene analysis database, expression pro les of the differentially expressed genes (DEGs) described at different stages of keloid development were utilized to explore the role of hub genes and immune-in ltrating cells in the development of this disease. To better describe the tumor gene expression characteristics of keloid development, we collected the tissue samples from keloid patients and divided them into three groups based on the keloid's development stage. We utilized a bioinformatics approach to identify potential tumor immune gene targets for keloids.

Patients and grouping
The clinical study protocol was reviewed and approved by the Medical Ethics Committee of Peking Union Medical College Hospital in Peking, China. All of the methods were performed in accordance with the relevant guidelines and regulations. Written informed consent was given by all patients. In all, 8 patients (5 keloid patients without in ammation other than the keloid mass and 3 keloid patients with in ammation besides the keloid mass) were randomly selected from the department of plastic surgery at Peking Union Medical College Hospital between January 2019 and March 2020. The information of keloid patients (gender, age) was collected, and the keloid's condition was assessed by modi ed Vancouver scar scale (mVSS) 31 according to the pigmentation, vascularity, pliability and height of the keloid lesions (Table 1). There were 4 male patients and 4 female patients in each group; their ages ranged from 18 to 37 years. All of the keloids were diagnosed and con rmed by an experienced plastic surgeon through pathological examination. No patients had any systemic disorders, were taking drugs, or were receiving other treatments that might affect the study results. All of the samples were taken from the chest region. Keloid samples were taken from the removed keloid tissue after operation. Normal skin samples were taken from normal skin tissue which had to be removed during keloid resection.
In ammatory tissue samples were taken from the in ammatory skin tissue after the acute stage of the three keloid patients with in ammatory lesions (Figure 1), which occurred in addition to the keloid mass and had to be removed at the same time during keloid resection.
The collected samples were divided into three groups: keloid tissue samples (K group), normal skin tissue samples from keloid patients (N group), and in ammation skin tissue samples after the acute stage from keloid patients who has in ammatory lesions besides keloids (I group). Table 1 The expression levels of tumor immune-related genes Oncomine Immune Response Research Assay kit (OIRRA, Thermo Fisher, USA) was used to analyze the tissue expression of tumor immune-related genes. The Oncomine Immune Response Research Assay kit is an RNA-based multi-target sequencing kit that can measure the expression levels of 395 genes related to immune response and tumor immunity in 36 functional annotation groups, consisting of genes related to lymphocyte regulation, cytokine signals, lymphocyte markers, checkpoint pathways, and tumor characteristics.
Intra-group data repeatability test The R programming language was used to provide software and operating environments for statistical analysis and graphic drawing. Principal component analysis (PCA) is a common method for sample clustering. It is usually used for gene expression, diversity analysis, and other sample clustering based on various information. Sample cluster analysis was used to test the repeatability of the data within the group.

DEGs identi ed and functional annotation
The statistical signi cance value was set as a p value < 0.05 and log2 fold change (FC) > 1.5. R package was used to create a volcano plot. Three data sets were utilized for lter DEGs. The Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 comprises a full knowledge base; this study utilized the sixth version of our original web-accessible program. DAVID [10] provides a comprehensive set of functional annotation tools for investigators to understand biological meaning behind an extensive list of genes. Also, it provides gene functional classi cation, gene ID conversion, and Gene Ontology (GO) term enrichment analysis to highlight the most relevant GO terms associated with a given gene list. The Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) is a database which help to understand the advanced functions and biological systems (such as cells, organisms, and ecosystems), generated from molecular-level information, especially large molecular datasets. The utility database resources were derived from genome sequencing and other highthroughput experimental techniques. These databases are often used to nd signal pathways for gene enrichment. GO standardized the description of gene products from functions, participating biological pathways, and localization in cells.
Construction and analysis of protein-protein interaction network, Hub gene analysis, and mining STRING [11] is a system that searches the known protein-protein interactions (PPIs) or predicts such interactions, which include both direct physical interactions and indirect functional correlations between proteins. After importing DEGs into the search tool, a PPI network is obtained. PPI analysis provides a new research basis for the pathophysiology of skin tissue. The minimum requirement for constructing the network is a medium credibility with an interaction score > 0.4. Cytoscape (version 3.6.1) is software that graphically displays networks, analyzing and editing them. Based on the topology principle, Cytoscape's plug-in molecular complex detection (MCODE) (version 1.5.1) detects densely connected regions of large PPI networks that may represent molecular complexes. The method is based on vertex weighting by local neighborhood density and outward traversal from a locally dense seed protein to isolate the dense regions according to given parameters. The algorithm has the advantage over other graph clustering methods of having a directed mode that allows ne-tuning of clusters of interest without considering the rest of the network, and allows examination of cluster interconnectivity, which is relevant for protein networks. Cytoscape software originally built a PPI network diagram. Second, MCODE [12] identi es the most important modules in the network map. In this study, criteria for MCODE analysis were degree cutoff = 2, MCODE score > 5, maximum depth = 100, k score = 2 and node score cutoff = 0.2. When the degree is set (degrees ≥ 10), the hub gene is excavated.

CIBERSORT
Determining the proportion of in ltrating immune cells during keloid development is an important aspect of studying the role of the keloid tumor immune microenvironment. The purpose of deconvolution is to try to infer the proportion information of immune cells from complex tissues. CIBERSORT [13] deconvolution algorithm is a machine learning method based on linear support vector regression and is highly robust to noise. The algorithm is superior to other methods in terms of noise, unknown mixture content, and closely related cell types. Through CIBERSORT deconvolution algorithm, the gene expression information is converted into immune cell proportion information, and through visual analysis of histograms, heat maps, and violin diagrams, one can isolate immune cells with signi cant differences.

Statistical analysis
Results are expressed as mean ± standard error of the mean. Statistical analysis was carried out using SPSS software version 21.0 (IBM corporation, Armonk, NY, USA), and p < 0.05 was deemed to indicate a statistically signi cant difference.

Evaluating data quality
To further validate the intra-group data repeatability, we used principal component analysis (PCA) to verify the repeatability of the data within the group. According to the PCA analysis, we found that the data in I group and N group were repeatable ( Figure 2A). PCA showed that the distance between the samples in the same group was very similar, and the distance between the samples in different compared groups was far in the principle component 1 (PC1) dimension (Figure 2A-2B). Due to the same patient tissue sample origination, there were partial data intersections between N group and K ( Figure 2C).

DEGs identi ed between the three groups
To identify the differences between N group, I group, and K group, we developed scatter plots (Figure 2A-2C) and volcano plots ( Figure 2D-2F). We identi ed 68 DEGs upregulated between I group and N group, and 41 DEGs were downregulated (Figure 2A and 2D). There were 46 DEGs upregulated between I group and K group, and 53 DEGs were downregulated ( Figure 2B and 2E). There were 21 DEGs upregulated between N group and K group, and 23 DEGs downregulated ( Figure 2C and 2F). Identifying the DEGs between the three groups was used to further identify hub genes.

Functional annotation of DEGs by GO and KEGG analyses.
A cutoff of log2 fold change > 1.5 and p < 0.05 in DEGs were the screening criteria used to perform GO and KEGG analysis. GO analysis revealed that there were many variations of BP between I group and N group, including regulation of lymphocyte activation and T-cell activation ( Figure 5D). The variations in DEG cell components (CC) were markedly enriched in on the external side of the plasma membrane and the plasma membrane receptor complex ( Figure 5E). Variations of MF were markedly enriched in cytokine receptor binding and cytokine receptor activity ( Figure 5F). GO analysis revealed that there were many variations of BP between I group and K group, including T-cell activation and chemokine response ( Figure  6A). The variations in DEG CCs were markedly enriched on the external side of the plasma membrane and early endosome ( Figure 6B). Variations of MF were markedly enriched in cytokine receptor binding and receptor ligand activity ( Figure 6C). GO analysis revealed that there were many variations of BP between N group and K group, including regulation of lymphocyte activation and leukocyte cell−cell adhesion ( Figure 6D). GO analysis also found that the variations in DEG CCs were markedly enriched on the external side of the plasma membrane and early endosome ( Figure 6E). The variations of MF in DEGs were markedly enriched in protease binding, integrin binding, and cytokine receptor binding ( Figure 6F). Analysis of the KEGG pathway between I group and N group revealed that all of the DEGs were primarily enriched in cytokine−cytokine receptor interaction, and viral protein interaction with cytokines and cytokine receptors ( Figure 6G). The result is the same as in I group and K group ( Figure 7A). Analysis of the KEGG pathway differences between N group and K group revealed that all of the DEGs were primarily enriched in rheumatoid arthritis and the intestinal immune network for IgA production ( Figure 7C).

The landscape of immune in ltration in keloid development
The molecular processes and the appearance of immune cells types were dynamic in different phases, which may be identi ed as a targeted therapeutic strategy. Our research aimed to describe differential tumor immune gene expression and immune cell composition in keloid development. To elaborate the complicated process of keloids, we collected three groups of tissues at the different stages of the disease. Owing to technical limitations, the landscape of immune in ltration in keloid development has not been entirely revealed before, especially for those with low abundance cell subpopulations. Using the CIBERSORT algorithm, we revealed the differences in immune in ltration between I group, N group, and K group of immune cell subpopulations.
The results obtained from I group and N group can re ect the tissue reaction process stimulated by in ammatory factors at the cellular level in keloid patients ( Figure 7C). Compared with tissue from I group, tissue from N group and K group generally contained a higher proportion of M2 macrophages ( Figure 7F), but the difference was not signi cant ( Figure 9A, p = 0.133; Figure 9B, p = 0.19). Figure 7D and Figure 7G summarized the difference between I group and K group. There was signi cant difference in M1 macrophages between these two groups ( Figure 9B, p = 0.03). The proportion of M1 in I group was higher than in K group. Figure 7E and Figure 7H summarized the difference between N group and K group. There were no signi cant subpopulation variations when taking the distinct group-bias clustering in PCA analysis into consideration between N group and K group ( Figure 2C). However, the proportion of CD8 + T cells varied signi cantly between N group and K group ( Figure 9C, p = 0.048). Furthermore, the proportions of different subpopulations of immune cells were moderately to strongly correlated ( Figure 8A-8C). Taken together, these results indicated that there was characteristic immune cell in ltration in different keloid evolution stage, which might have important pathological meaning in keloid development.

Discussion
Keloids are a dis guring broproliferation disorder that dramatically impairs the quality of life of the affected individual. Even though there have been many studies about keloids, they remain the most challenging skin condition to treat. The theory of tumor immune microenvironment states that the occurrence and development of tumors are the result of the joint action of various cells and molecules in the microenvironment of the tissues where tumor cells are located [14]. Keloids have certain tumor characteristics. We speculate that the imbalance of skin immune microenvironment caused by genetic, immune, in ammatory, and other factors may be key to keloid development. Among them, abnormal expression of immune-related genes is the root cause, and immune cell regulation plays a key role in the development of keloids. These factors result in keloid tumor-like scar hyperplasia formation. Aberrant gene expression often originates from congenital risk, triggered by stimulating factors. Expression pro les of mRNAs, lncRNAs, and circRNAs were altered in keloid tissue, which may partly contribute to the etiology of keloids by affecting several signaling pathways related to keloid development [15]. DEGs can be used as biomarkers for diagnosis or as therapeutic targets for keloids. Some biomarkers are upregulated in keloids, including broblast-speci c protein 1 (FSP1), vimentin, TGF-β1, and Smad3 phosphorylation, which were noted in keloid tissues [16]. Expression of miRNAs in keloid tissue was shown to be signi cantly different from that of normal tissue. Among the miRNAs involved in keloid pathogenesis, miRNA-21, miRNA-141-3p, miRNA-181a, and miRNA-205 were thought to upregulate broblast proliferation and depress apoptosis of keloid-derived broblasts through the PI3K/Akt/mammalian target of rapamycin (mTOR) signaling pathway [17]. Drugs that indirectly target the biochemical microenvironments of keloids include growth factors (e.g., TGF-bs, bFGF, and VEGF), immunomodulators (e.g., tacrolimus, trehalose, IFN, and imiquimod), and anti-in ammatory drugs (e.g., IL-10 and IL-6) [18][19][20][21][22][23][24]. Based on the above research, we believe that keloid patients have abnormal gene expression compared to ordinary people.
To better understand keloids pathogenesis, we collected three groups of tissues from patients in different disease stages to describe the development of keloids. To identify abnormal gene expression, we collected normal skin tissue from keloid patients, which we described as N group. To determine the gene expression changes, we collected in ammatory lesion tissue from keloid patients, which we described as I group. Based on the tumor immunological gene study, we identi ed hub genes from N group to I group. These genes included CCR1, CCR2, SELL, IL10, CCR7, CD40LG, CD69, CXCL8, IL6, and CXCL9. Current hypotheses of pathogenesis classify keloids as an entity of aberrant brosis and tumor characters.
Hyperactivation of the MCP-1/CCR2 axis reportedly causes brosis in the liver, cirrhosis, atherosclerosis, and lung brosis [25][26][27]. Expressions of MCP-1 and its receptor CCR2 in keloid lesions were increased in the keloid tissues [25]. IL-10 was shown to be able to signi cantly inhibit the proliferation of keloid broblasts [28]. Moreover, proin ammatory factors, such as interleukin (IL)-1α, IL-1Ra, IL-1β, IL-6, and tumor necrosis factor-α are upregulated in keloid tissues, which suggests that, in patients with keloids, proin ammatory genes in the skin are sensitive to trauma [1,29]. Hub genes between I group and K group included IL10, ITGAM, ITGAX, IL2, IL4, IL6, IL13, IL17A, FOXP3, and CD86. The roles of several genes in keloid development, including CCR1, SELL, CCR7, CD40LG, CD69, CXCL8, ITGAM, ITGAX, CD86, and CXCL9, need to be veri ed. These genes may play an important role in keloid formation after stimulation by in ammatory factors. We need to further validate these hub genes to reveal the molecular mechanism of keloid development. This may help to nd a new therapeutic target.
The occurrence and development of keloids involve complex molecular mechanisms. We used GO and KEGG analysis to better understand the mechanisms of keloids. A study by Dohi T et al. identi ed that the soft skin surrounding keloids is exposed to high mechanical strain that correlates with increased expression of the caveolin-1/rho signaling via the rho kinase mechanotransduction pathway, which may lead to keloid progression. A study by Huang H et al. identi ed that upregulated mRNAs were involved in cell proliferation, cell growth, and tissue repair, and downregulated mRNAs were involved in apoptosis [30]. In our study, GO analysis revealed that there were many variations of BP between I group and N group, including regulation of lymphocyte activation and T-cell activation. The result is similar to that found in I group and K group, which may play an important role in the initiation and formation of keloids.
Variations of MF were markedly enriched in cytokine receptor binding and receptor ligand activity, which are essential for signaling pathway transmission.
In this study, analysis of the KEGG pathway between I group and N group revealed that all of the DEGs were primarily enriched in cytokine−cytokine receptor interaction and viral protein interaction with cytokine and cytokine receptor, re ecting that these pathways may affect keloid formation. A study by Zhong L et al. identi ed that target genes were associated with the MAPK signaling pathway and HIF-1 signaling pathway [31]. Our study provides a new eld for the molecular study of keloids.
The mechanistic details of keloid formation remain poorly understood. Given that the immune system is engaged in skin lesion repair, we explored cell-targeted analysis in this study. We utilized CIBERSORT to transform gene expression pro les to cell composition of complex tissues. There are immune cells, broblasts, endothelial cells, and an abundant collection of cytokines, chemokines, and growth factors in keloid tissue. These components and their complex interactions form a tumor-associated microenvironment, and tumor cells may use immune cells to seek bene ts of growth, invasion and transfer. The main host cells recruited and activated in the tumor microenvironment are various immune cells. These immune cells have important prognostic relevance due to the dual role of the immune system in promoting and suppressing tumors. Searching for immune cells that may play a major role in keloid development will be bene cial to clinical application of cell therapy as a target. M1/M2 proportion describes the two major and opposing activities of macrophages [32]. M1 macrophages secrete proin ammatory factors, chemokines, and presenting antigens full-time, participate in the positive immune response, and play a role in immune surveillance [33]. M2 macrophages are an important immune cell in downregulating the immune response [33]. Remarkably, the molecules primarily responsible for these "Fight" (NO) or "Fix" (Ornithine) activities both arise from arginine, and via enzymatic pathways (iNOS and arginase) that down regulate each other [32]. The names M1 and M2 were chosen because M1 and M2 macrophages promote T1 and T2 responses, respectively. Products of T1 and T2 responses (e.g., IFNγ and IL-4) also downregulate M2 and M1 activity, respectively [34]. Thus, M1/M2 proportion demonstrated the importance of innate immunity, and how it is linked to adaptive immunity in a beautifully counterbalanced system [35,36].
M2 macrophages have only a weak antigen presentation function and play an important role in immune regulation by secreting congruent cytokines, such as IL-10 and TGF-β. They are an important immune cell in downregulating immune response [33]. In this study, we identi ed a higher proportion of M2 macrophages in N group than in I group. M2 macrophages may lead to immunosuppressive regulation. Its high proportion in skin tissue in keloid patients may indicate the abnormal gene expression in the skin, which leads to a state of immunosuppression. This also means that the skin of keloid patients is more vulnerable to infections that cause in ammation and keloid formation. However, this difference was not signi cant probably due to insu cient sample size and still need to be veri ed by comparing with the skin from normal person.
Our research showed that there was a signi cant difference in M1 macrophages between I group and K group. The proportion of M1 in I group was higher than in K group. It is supposed that the state of immunosuppression could be activated by in ammatory factors. M1 macrophages play an important role in activating the immune state in keloid patients during infection and in ammation. This may indicate the presence of immune activation in the in ammation stage of keloid development. A study by Jin Q et al. identi ed that macrophages in keloid tissues were polarized toward the M2 subtype [8]. In our study, the proportion of M2 in K group was higher than in I group, although the difference was not obvious. Keloid tissues presented signi cantly elevated in ltration by CD14+ macrophages [25]. Myo broblast proliferation and heterogeneity are supported by M2 macrophages during skin repair [37].
This means that M1/M2 proportion in different disease stages may play a major role in keloid development.
In addition, the proportions of CD8 + T cells varies signi cantly between both N group and K group.
Contrary to our ndings, Jin Q et al. identi ed that keloid tissues presented higher in ltration by CD3 + T cells, of which the majority were CD4 + T cells. This difference may be caused by different tissue sample origination, and needs to be clari ed with large sample analysis.
Despite the rigorous bioinformatics analysis in this study, there are still some limitations. First, the sample size of our study was small, which might result in some deviations in the results. Second, we only conducted bioinformatics mining and lack experimental veri cation, which will be carried out in the future.

Conclusion
Taken together, we traced multiple hub genes involved in keloids formation and immune cell in ltration character during keloid development. Several genes may play an important role in the development of keloids that have not been studied in literature before, including CCR1, SELL, CCR7, CD40LG, CD69, CXCL8, ITGAM, ITGAX, CD86, and CXCL9. BP, CC, and MF that associated with tumor immune gene expression in different keloid development stages were described in detail using bioinformatics. Immune in ltration in different disease stages may play an important role in keloid development. M1/M2 macrophage proportion and CD8+T cells may be the key immune cells in the keloid development process.

Consent for publication
Not only did all of the patients participating in the experiment sign the hospital's surgical informed consent and photographic informed consent, but they also signed the informed consent for the use of tissue samples for donation under voluntary principles.

Availability of data and materials
The datasets during and analyzed during the current study available from the corresponding author on reasonable request.

Competing interests
No competing interests are involved in this research.

Funding
This study was supported by The National Natural Science Foundation of China (81871538).  Figure 1 In ammatory tissue samples were taken from the in ammatory skin tissue after the acute stage of the three keloid patients with in ammatory lesions green arrow , which occurred in addition to the keloid mass and had to be removed at the same time during keloid resection.   The hub genes were identi ed from the PPI network in N group and K group. functions (MF) of DEGs between I group and N group. Cutoff for log2 fold change > 1.5 and p < 0.05 as screening criteria.