Underlying Genes and Molecular Mechanism of Keloids Investigated by Integrated Bioinformatics Analysis

Background: We aimed to identify the overlapping differentially expressed genes (DEGs) of keloids distinguished from normal scar and normal skin and relevant underlying mechanism using integrated bioinformatics methods. Methods: The expression proles of 18 keloid samples, 7 normal skin and 5 normal scar, were obtained from the GSE7890, GSE44270, GSE92566, and GSE3189 datasets in the Gene Expression Omnibus database. DEGs were identied using the LIMMA package in R. Gene ontology (GO) functional enrichment analysis was performed using the R software. A DEG-associated protein–protein interaction (PPI) network was constructed using STRING and MCODE was used for module analysis of the PPI network. Moreover, the hub genes were veried by qRT-PCR. The predicted DEGs, their regulatory miRNA and TF regulation network was analyzed using miRnet. Results: A total of 978 common DEGs were identied in the keloid samples. Genes with more than 45 interaction degrees, including neuropeptide Y (NPY), opioid receptor mu 1 (OPRM1), cholinergic receptor muscarinic 2 (CHRM2), and proopiomelanocortin (POMC), were found in the PPI network. Hsa-miR-335 and Sp1 as upstream-regulators regulated CHRM2, NPY, and POMC. Functional enrichment analysis revealed that hub genes were commonly enriched in the “G protein-coupled receptor signaling pathway” GO_BP term Conclusion: Taken together, CHRM2, NPY, POMC, and OPRM1 potentially have crucial roles in keloid disease. Furthermore, miR-335 and Sp1 are potential targets for preventing keloid formation.


Data Resources
The mRNA expression pro les of GSE7890, GSE44270, GSE92566, GSE3189 (all homo sapiens) were acquired from Gene Expression Omnibus (GEO) database. Five normal scars and ve keloid broblast samples were extracted from GSE7890. Three normal skin and nine keloid broblast samples were downloaded from GSE44270. Four keloid samples were derived from GSE92566. Other four normal skin samples were derived from GSE3189. Samples of both GSE7890 and GSE92566 datasets were located on GPL6244 (Affymetrix Human Gene 1.0 ST Array). Others belong to GPL570(Affymetrix Human Genome U133 Plus 2.0 Array) and GPL96(Affymetrix Human Genome U133A Array), respectively.

Data preparation and screening of DEGs
Downloading series matrix le and platform le for conversion into gene symbol from NCBI. Using R project for statistical computing (Version 4.0.0) that is inserted by both LIMMA (https://www.bioconductor.org/packages/release/bioc/html/limma.html) package to run a series of procedure, as followed ID conversion, merging, batch analysis for normalization and expression calculation.
The nal calculation of genes which were matched with multiple probes was the average of those probes value.
After mentioned operation, the screening of DEGs was carried out via LIMMA package in R. Three groups, including A [keloid samples(n=18) from excision surgery], B [normal scar samples(n=5) from adults], C [normal skin samples(n=7) from elective plastic surgery and melanoma patients], were de ned. To identify much more potential mRNA that are related with the pathogenesis of keloid, keloid groups(A) were compared with B, C groups, respectively. Then, we presented the common DEGs in the online Venn diagrams tool (http://bioinfogp.cnb.csic.es/tools/venny/index.html and http://bioinformatics.psb.ugent.be/webtools/Venn/), which were for subsequent analysis. The DEGs were obtained base on the lter condition of p-value 0.05 and |log2 fold change (FC)| 1, which were depicted by volcano plot.

Functional Enrichment Analysis
The Gene Ontology (GO) functional enrichment analysis of DEGs were accomplished by the cluster Pro ler R package.
For detailed information of DEGs, we classi ed screening genes by biological process (BP), molecular function (MF), cellular component (GC) and the corresponding signaling pathway with the GO (gene ontology) analysis [11] and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis [12], which displayed in ClueGO (Version2.5.5)plugin of Cytoscape (Version 3.7.2) [13].

PPI Network Analysis and Subnet Module Analysis
Applying with the o cial website (http://string-db.org/) of Search tool for the retrieval of interacting genes/proteins (STRING) acquired the interaction network among the DEGs-encoded proteins and the medium con dence was required for ≥0.9 [14]. The most signi cant PPIs were presented in Cytoscape software (version: 3.7.2, http://cytoscape.org/) [15] with MCODE application and acquired the most signi cant module under k-core>10 circumstance. The genes located in the same module tend to generate proteins that possess the same or analogous function, moreover, these proteins integrate into a module in a same biological role. Usage of the MCODE algorithm to calculate the score of each module. With the higher score of the module showed the more closed interactions and enrichments. For obtaining hub genes, PPI network were analyzed using R.
Prediction of hub-genes with miRNA and TF mutual regulating networks.
The miRnet, a visual network-based analysis, which could accomplish relationship of miRNA and TF with targets according to the microRNA target base(miRTarBase) and microRNA records(miRecords). After establishing regulation network, cytohubba was used to select the highest degree of miRNA and TF.

Quantitative real time polymerase chain reaction(RT-qPCR)
Total RNA was isolated and performed reverse transcription to acquire stable cDNA that was used to subsequent experiment. The expression of NPY, OPRM1, CHRM2, POMC was quanti ed by the 2( -▲▲Ct ) method. Results:

Selection of DEGs
The DEGs from two sets were screened at p<0.05 and | log FC | >1 standard (Fig.1A, B). A total of 2714 genes were up-regulated and 926 were down-regulated genes as displayed in the volcano plot. The common genes screened from the two sets included 952 up-regulated genes and 26 down-regulated genes as displayed by the Venn tool (Fig.1C).
All the genes of mcode3 and mcode4 participated in "extracellular space" (GO:0005615) and "intracellular part" of the GO_CC terms, respectively.
Identi cation of miRNA-target and miRNA-TF interaction network Furthermore, we detected the possible correlated miRNA and TF through online analysis tool-miRnet for DEGs.

Discussion:
In this study, the common 978 DEGs(952 up-regulated and 26 down-regulated) in keloid samples were markedly differentiated from normal skin and normal scar samples. According to the analysis of the PPI network, the core genes were as followed: GNG3, GNG13, POMC, ADCY8, NMUR1, SST, CCR3, GCG, CHRM2, NPY, OPRM1, PDYN and PYY. Meanwhile, the miR-335-5p and Sp1 were identi ed as the main upstream factors, of which Sp1 was found to be associated with NPY, OPRM1 and miR-335-5p modulated three central genes that were NPY, POMC, CHRM2. Except for PYY gene, these core genes were commonly enriched in G protein-coupled receptor signaling pathway" in GO_BP terms. It is noticed that NPY, POMC, CHRM2, OPRM1 all belong to the mocode1, which underlines the similarities of their functions.
In the keloid microenvironment, the mechano-physiological conditions are important, they are involved in mechanical tension and in ammatory tension and a vicious circle was formed in between. The mechanical information needs to be converted into a biological signal through the cell membrane mechanoreceptors.
Besides integrins, there are non-integrins for accomplishing the transformation process such as G-protein coupled receptors [16]. In the present study, both OPRM1 and CHRM2 belong to a superfamily of the G-proteincoupled receptors, inferring their latent correlation with the mechanical tension. Moreover, these receptors were also involved in growth and progression of tumor [17,18]. In the view of the abnormal proliferation of keloid broblasts is the major reason for keloid formation, the similar role of receptor genes is speculated in keloid.
Altered balance between proangiogenic and antiangiogenic is required for tumor growth beyond a certain size. Vascular endothelial growth factor (VEGF), periostin, and endostatin have been shown to be related to new vessel formation in keloid lesions [19,20]. We found that NPY, which encodes another direct angiogenic stimulator, was up-regulated in the keloid samples in our study. NPY-stimulated VEGF production and secretion was found to contribute highly to angiogenesis activity in human breast cancer [21] and NPY also was identi ed as a promoter of prostate and breast cancer, affecting the proliferation and migration of cells [22,23]. Furthermore, NPY was shown to participate in brogenesis of hepatic stellated cells, thereby contributing to hepatic cancer [24].
The deregulation of extracellular matrix deposition is one of the pathological processes in keloids. POMC mRNA was detected in keloid-derived broblast, which is in agreement with our results [25]. As the precursor of various active peptides, POMC-derived products were found to have biological roles under cytokine stimulation in the regulation of extracellular matrix deposition and in ammation [26],suggesting their potential functions during the development of keloids.
Among the 2049 predicted miRNAs, miR-335 had the highest degrees and targeted DEGs in the regulation network. MiR-335 has been veri ed as both a tumor suppressor and tumor promoter in various cancers [27].
One of the inhibition mechanisms of miR-335 is the activation of tumor suppressor p53 by alerting Rb1 to repress cell proliferation [28], and the antagonistic effect of miR-335 on miR-21 was shown to be mainly a protumorigenic mechanism [29]. miR-21 involvement in keloids has also been reported [30], suggesting that miR-335 is involved in the underlying mechanism and neuropeptides, including NPY, POMC, and CHRM2, are involved downstream. This axis mechanism requires further validation. SP1 is a well-known TF that is involved in keloid pathogenesis mainly by regulating the extracellular matrix process of keloids [31]. However, a more detailed understanding of the underlying mechanism is required. We found that Sp1 formed an interaction network with NPY and OPRM1 in the TF-DEG network.

Conclusion:
Taken in combination,the common 978 DEGs (952 up-regulated DEGs and 26 down-regulated DEGs) were con rmed by comparing keloid samples with normal skins samples and normal scars samples. As for screened core genes from PPI network, especially for NPY, POMC, CHRM2, OPRM1 that was supported by RT-qPCR experimental validation. Aimed at four genes, we established a interaction network to identify the upregulators including miRNA and TFs, moreover, GO functional enrichment analysis provided clues to speculated the "G protein-coupled receptor signaling pathway" GO_BP is highly related with the molecular of keloids development. Furthermore, certain observation should be made by experimental validation.