Identification and Validation of Novel Diagnostic Biomarkers for Keloid Based on GEO Database.

INTRODUCTION: Keloid is a pathological scar type, which invades normal surrounding tissue without self-limiting to cause pain, itching, cosmetic disfigurement, etc. Knowledge of the molecular mechanisms underlying keloid remains unclear. This dilemma leads to no biomarker available for diagnosis. Thus, to seek accurate diagnosis, biomarkers are necessary for keloid diagnosis to help control its incidence.
METHODS: Gene Expression Omnibus (GEO) database was used to select differentially expressed miRNAs (DE-miRNAs) in GSE113620. miRTarBase miRNA-target tools were used to predict the interactions between miRNAs and their target mRNAs. Target mRNAs that were differentially expressed in keloid were selected by analyzing differentially expressed genes (DEGs) in GSE44270 and GSE92566. PPI network analysis, gene enrichment analysis, cell-specific and tissue-specific expression analyses of DE-target mRNAs were conducted. RT-PCR analysis was conducted to validate our results.
RESULTS: Three novel miRNAs (miR-30b-5p, miR-212-3p, miR-149-5p) and five target mRNAs (SIX1, CCNA2, CCNB1, FOXM1, RUNX2) were identified as potential biomarkers for keloid patients. Additionally, the potential functions of those miRNAs-mRNAs pathways were analyzed. DISCUSSION: These findings of keloid-related miRNAs, mRNAs, and miRNA-mRNAs regulatory networks may provide insights into the underlying pathogenesis of keloid and serve as potential biomarkers for keloid diagnosis.

Introduction

Keloid is excessive fibrosis of the skin that extends beyond the area of injury and does not regress.1 Keloid can occur in the joints and mouth after several years of severe injury, including burns, chemical injury, wound, and surgical incision.2 Keloids on the joints affect the quality of life, physical, mental, and social health.3 Knowledge of the molecular mechanisms underlying keloid remains unclear. This dilemma leads to no biomarker available for diagnosis, and, therefore, result in the diagnosis of keloid in clinical practice is highly dependent on the experience of the clinician. Thus, to seek an effective biomarker may offer an alternative to the current clinical pain point.

MicroRNAs (miRNAs) are evolutionarily conserved short (20–23 nucleotides) noncoding RNAs with no protein-coding potential but can negatively or positively regulate genes.4 Previous studies have shown that miRNAs can inhibit gene expression by binding to the 3‘-untranslated regions (UTR) of target mRNAs, thus inhibiting mRNA translation or promoting mRNA degradation.5,6 Several studies have revealed that miRNAs are involved in various biological pathways, such as cell proliferation, apoptosis, cell cycle, differentiation, metabolism, and abnormal expression of miRNAs, which may impair a variety of cell functions, leading to cancer and other diseases.7 Several studies have shown that miRNAs are biomarkers for keloid.8 Therefore, it is important to reveal the biological function of miRNAs in keloids to analyze the mechanism of fibroblasts in keloids effectively. This study aimed to investigate several potential biomarkers for distinguishing keloid patients from non-keloid patients.

Gene Expression Omnibus (GEO) database was used to identify DE-miRNAs and their target mRNAs differentially expressed in keloid and control groups. A DE-miRNA and differentially expressed target mRNAs (DE-target mRNAs) network was then constructed after a comprehensive analysis. Further, protein–protein interaction (PPI) network analyses, functional enrichment analysis, tissue-specific and cell-specific gene expression analysis were performed. Validation analysis was conducted using RT-PCR experiment. Three potential miRNAs and five target mRNAs were identified as significant predictors of keloid. This study may provide novel diagnostic biomarkers for keloid.

Materials and Methods

Data Source and Identification of DE-miRNAs/DE-mRNAs

MicroRNA expression profile data of GSE113620 were obtained from Gene Expression Omnibus (GEO; ) database, which is based on GPL19117 [miRNA-4] Affymetrix Multispecies miRNA-4 Array. RNA-seq detection and miRNA expression data were analyzed in serial skin biopsies from the same site (baseline and 6 weeks post-injury) of keloid-prone (n=8) and healthy matched control individuals (n=6). Differential expression analysis was conducted in four comparison groups: 1) Keloid-prone subjects vs healthy subjects in baseline, 2) Keloid-prone subjects vs healthy subjects 6 weeks after wounding, 3) Keloid-prone subjects in baseline vs keloid-prone subjects 6 weeks after wounding, and 4) Healthy subjects in baseline vs healthy subjects 6 weeks after wounding. The mRNA expression profile data of GSE44270 and GSE92566 were obtained from the GEO database (which are based on GPL6244 and GPL570, respectively) to identify common mRNAs between the target mRNA of DE-miRNA and DEGs in GSE44270 and GSE92566. For GSE44270 analysis, keloid scars were excised from patients undergoing scar excision surgery. Normal skin samples were isolated from patients undergoing elective plastic surgery. This study had six experimental groups, including nine cases of keloid fibroblast (KF), nine cases of keloid keratinocyte (KK), four cases of keloid non-lesional fibroblast (NKF), four cases of keloid non-lesional keratinocyte (NKK), three cases of normal fibroblast (CF) and three cases of normal keratinocyte (CK). For GSE92566 analysis, transcriptional profiling was conducted using Affymetrix HGU133 2.0 plus arrays for three groups, including large chronic keloids (n=3), adjacent non-lesional keloid samples (n=3), and a newly formed keloid lesion. The Affymetrix CEL files and the probe annotation files were obtained from the GEO database. Meanwhile, all the samples were preprocessed with background correction. The Bioconductor package () (R studio Version 1.1.383) was used to select DE-miRNAs. The |log2FC (fold change) |>1 and FDR <0.05 were used as the cut-off criteria for selecting DE-miRNAs. “ggplot” package was used for data visualization.

Identification of Target mRNA of DE-miRNAs

Identifying target mRNAs of miRNAs is a necessary step for analyzing the regulatory mechanism and function of miRNAs. Herein, miTarbase was used to identify the target mRNAs of 61 DE-miRNAs. The DE-miRNAs, with more than five target genes (support type: functional MTI), were selected. Cytoscape software was then used to visualize the regulatory networks of the miRNA–mRNA pairs.9 The common mRNAs between the target mRNAs and the DEGs in GSE44270 and GSE92566 were then identified. An online Venn tool () was used for Venn diagram analysis.10

Protein–Protein Interaction (PPI) Network Analysis

A String database () (a database of known and predicted protein–protein interactions) was used to construct a PPI network to analyze the interactions between target mRNA and peripheral genes. The nodes and lines in the network diagram represent target genes and their interactions, respectively, and the color of the lines represents the strength of the interaction.

Gene Ontology (GO) Annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analyses of the DEGs

The Database for Annotation, Visualization, and Integrated Discovery (DAVID) was used for GO annotation and KEGG pathway enrichment analyses of the selected mRNA cluster. GO annotation analysis included biological processes (BPs), cellular components (CCs), molecular functions (MFs). R studio “ggplot” package was used for data visualization.

Tissue-Specific Gene Expression Analysis

An online resource BioGPS (, was used to analyze the tissue-specific expression of the target mRNAs. Herein, a gene was defined as highly tissue-specific if the level of tissue-specific expression was 10 times that of the median. The top three tissues with the highest expression levels were recorded if a gene was expressed more than 10 times in multiple tissues.

Bgee is a database that compares gene expression patterns across a wide range of animal species, generated by integrating multiple data types (RNA-SEq, Affymetrix, in situ hybridization, and EST data). It is based entirely on selected healthy wild-type expression data (eg, no knockout, no treatment, no disease) to provide a comparable reference for normal gene expression. Thus, Bgee can detect the expression conditions of any single gene, adapted to any data type and species. Bgee version 14.1 includes 29 animals and is available at and the development pack BgeeDB.

Cell-Specific Gene Expression Analysis

HPA () is an online database that can be used to detect the expression of target genes in different cells, tissues and organs. The database is based on RNA sequencing (RNA-seq) analysis and immunohistochemical (IHC) data from previous studies.11 Herein, the HPA database was used to analyze the expression of target mRNAs of DE-miRNAs in multiple single cells of the skin. The single-cell RNA sequencing data in the HPA database were obtained from publicly available whole-genome expression data. The site contains all protein-coding genes in 12 cell-type groups and 192 cell-type clusters.11 It enables the user to study the expression of genes in different single-cell types and their expression in different single-cell types. The genes expressed in each cell type can be explored in interactive UMAP plots and bar charts, with links to corresponding immunohistochemical staining in human tissues.11

Sample Collection and Storage

Patients with keloids and healthy individuals were recruited from the Clinic of Department of Plastic and Reconstructive Surgery, Shanghai Ninth People’s Hospital. Keloids were diagnosed as previously described.12 Scar tissue from keloid patients and normal skin tissue from the inguinal area of healthy individuals were used for validation tests. Each collected sample was labeled and stored in liquid nitrogen until needed. The basic information of keloid patients and healthy donors is described in Table 1.

Table 1

General Information of Patients

	Keloid Patients	Normal Control
Number	16	20
Males/females (n)	7/9	4/16
Age (y), mean ± SD	29.3±11.4	29.62±15.22
Surgery/injury history	12	15
History of infectious diseases (HIV, viral hepatitis, etc)	2	4
Hormone therapy history	13	2
Laser treatment history	13	2
SCAR score (n)	11 (4) 10 (2) 9 (4) 8 (6)	N/A
Site of keloid	Lower trunk (1); upper trunk (14); face (1)	N/A

Total RNA Extraction and Quantitative RT-PCR

Total RNA was isolated from keloid sample and normal skin sample using TRIzol reagent (Invitrogen). Next, cDNA was synthesized using the PrimeScript RT reagent Kit (Takara Bio Company, Shiga, Japan). After reverse transcription of total RNA, quantitative RT-PCR was conducted to determine the expression of MAIT cells using SYBR Premix Ex Taq (TaKaRa, Dalian, China) on a Light Cycler 480 SYBR Green I Master (Roche Diagnostics, Basel, Switzerland). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) served as the control gene. The 2-ΔΔCt method was used to determine the gene expression levels. Primer sequences are shown in Table 2. All of the experiments were performed three times in triplicate.

Table 2

Information of Primer Sequences Used for RT-PCR

Gene Symbol	Primer Sequences
miR-30b-5p	Forward: 5′-GGGCTGTAAACATCCTACAC-3′
	Reverse: 5′-GGGCTGTAAACATCCTACA-3
miR-212-3p	Forward: 5ʹ-AGCATCCACGAGCAAGAGAC-3’
	Reverse: 5ʹ-GATGCTACTAGTGTGGCGGG-3’
miR-149-5p	Forward: 5ʹ-CCCTCATTCTGTGCCACACTCCAGCTGGG-3′
	Reverse: 5ʹ-TGGTGTCGTGGAGTCG-3′
U6	Forward: 5ʹ-GCTTCGGCAGCACATATACTAAAAT-3’
	Reverse: 5ʹ-CGCTTCACGAATTTGCGTGTCAT-3’
SIX1	Forward: 5ʹ-AAGGAGAAGTCGAGGGGTGT-3′
	Reverse: 5ʹ-TGCTTGTTGGAGGAGGAGTT-3′
CCNA2	Forward: 5ʹ-GGATGGTAGTTTTGAGTCACCAC-3′
	Reverse: 5ʹ-CACGAGGATAGCTCTCATACTGT-3′
FOXM1	Forward: 5′-TGCCCAGCAGTCTCTTACCT-3′
	Reverse: 5′-CTACCCACCTTCTGGCAGTC-3′
CCNB1	Forward: 5′- GCACTTTCCTCCTTCTCA-3′
	Reverse: 5′- CGATGTGGCATACTTGTT-3′
RUNX2	Forward: 5′- CACTGGCGGTGCAACAAGA-3′
	Reverse: 5′- ATGACGGTAACCACAGTCCCATC-3′
GAPDH	Forward: 5ʹ-GGAGCGAGATCCCTCCAAAAT-3′
	Reverse: 5ʹ-GGCTGTTGTCATACTTCTCATGG-3′

Statistical Analysis

All the data were analyzed using GRAPHPAD PRISM (version 6, GraphPad, La Jolla, CA, USA). Differences between groups were analyzed using Student’s t-test. The significance was defined at P < 0.001 (***), P < 0.01 (**), and P < 0.05 (*).

Results

Identification of DE-miRNAs

MicroRNA expression profile data of GSE113620 were obtained from the GEO database. A total of 164 DE-miRNAs (91 downregulated and 73 upregulated genes) were found in keloid-prone subjects in baseline vs keloid-prone subjects 42 days after wounding (Figure 1A). A total of 114 DE-miRNAs (38 upregulated and 76 downregulated genes) were found in healthy subjects in baseline vs healthy subjects 42 days after wounding (Figure 1B). Venn diagrams were developed to analyze common genes in different subsets, of which 69 genes were expressed differently only in the keloid-prone subjects in baseline vs keloid-prone subjects 42 days after wounding (Figure 1C). Those 69 DE-miRNAs were selected for subsequent analysis.

Figure 1

(A) Identification of DE-miRNAs. A total of 164 DE-miRNAs, including 91 downregulated and 73 upregulated genes, were found in keloid-prone subjects in baseline vs keloid-prone subjects 6 weeks after wounding. (B) However, 115 DE-miRNAs, including 38 upregulated and 76 downregulated genes, were found in healthy subjects in baseline vs healthy subjects 6 weeks after wounding. (C) Venn diagrams were developed to analyze common genes in different subsets, of which 69 genes were expressed differently only in the Keloid 0-Day vs Keloid 42-Day group.

Identification of Target mRNAs of DE-miRNAs

The target genes of the above 69 genes were identified in Mirtarbase, the target gene prediction database for miRNAs. Six DE-miRNAs, with more than five target genes (support type: functional MTI), were used to build a network of DE-miRNAs-mRNAs interactions using Cytoscape (Figure 2A). The six DE-miRNAs and their target mRNAs are shown in Table 3 (159 genes are listed, eight target genes overlapped). The 151 target mRNAs were then compared with the DEGs in two microarrays (GSE44270 and GSE92566) to identify the target mRNAs related to keloid. There were six comparison groups in GSE44270 and one in GSE92566 (KK vs NLK, KF vs CF, KK vs CK, KF vs NLF, NLK vs CK, NLF vs CF in GSE44270, and Keloid vs Control in GSE92566). There were no common DEGs in the KF vs NLF group, which is not shown in the figure Finally, among the 151 target mRNAs of DE-miRNAs, we found seven mRNAs that were differentially expressed in keloid by analyzing DEGs in GSE44270 and GSE92566. LIF and SIX1 were identified in KF vs CF; CCNA2, CCNB1, and FOXM1 were identified in KK vs NLK; DPYSL2 was identified in KK vs CK; RUNX2 was identified in Keloid vs Control (GSE92566). The final results were expressed as a Venn diagram (Figure 2B). The PPIs were mapped using the logical data from the STRING database to further assess the correlation between the above seven genes (Figure 2C). The network consisted of seven nodes and five edges, and the average local clustering coefficient was 0.524. The PPI network showed that SIX1, CCNA2, CCNB1, FOXM1, and RUNX2 formed a gene cluster regulated by three DE-miRNAs (Figure 2D). Those five genes were used as DE-target mRNAs in the subsequent analysis.

Figure 2

(A) MiRNA–target gene interactions. By searching the miRNA target gene prediction database, we found the target genes of DE-miRNAs. The predictions were verified by mirtarbase. Eventually, we selected six DE-miRNAs, whose target genes (support type: functional MTI) are more than five and build the network of miRNA–mRNA interactions using Cytoscape. Red node indicates the expression of gene is increased in keloid, while green indicates the expression of gene is decreased in keloid. (B) After selecting mRNAs of DE-miRNAs, we compared the above target mRNAs with DEGS in GSE44270 and GSE92566, and selected the overlapping mRNAs for future studies. Finally, we selected seven mRNAs that were DEGs in GSE44270 or GSE92566. (C) We plotted the PPI network diagram of the above seven genes to describe the connections among the seven genes. Moreover, LIF and DPYSL2 genes had no association with the other five genes, hence were abandoned in further study. (D) The connection of three DE-miRNAs and five DE-target mRNAs. Red node indicates the expression of gene is increased in keloid compared with control groups, while green node indicates the expression of gene is decreased in keloid compared with control groups.

Table 3

Target Genes of Identified DE-miRNAs

The miRNA–mRNA Network
Gene Symbol	Experiments	References (PMID)
Target genes of hsa-miR-181c-5p (Up regulated, log2FC= 1.671, adj P value=0.01)
NLK	Luciferase reporter assay	19585654
GATA6	Luciferase reporter assay	19585654
CDX2	Luciferase reporter assay	19585654
KRAS	Luciferase reporter assay/qRT-PCR/Western blot	20080834
NOTCH4	Luciferase reporter assay/qRT-PCR/Western blot	20080834
BCL2	Luciferase reporter assay/qRT-PCR/Western blot/	20162574
TRIM2	In situ hybridization/Luciferase reporter assay	21720722
SIRT1	In situ hybridization/Luciferase reporter assay	21720722
BTBD3	In situ hybridization/Luciferase reporter assay	21720722
BMPR2	Luciferase reporter assay	23352489
IL2	Reporter assay; qRT-PCR	21112091
DPYSL2	qRT-PCR/Western blot	23667495
MECP2	qRT-PCR/Western blot	23667495
RAP1B	Immunofluorescence/Luciferase reporter assay/qRT-PCR/Western blot	24573637
PRKCD	Luciferase reporter assay	26629004
PRAP1	qRT-PCR/ELISA/Luciferase reporter assay/Western blot	25554417
MST1	Immunoprecipitation/Luciferase reporter assay/Western blot	26561204
LATS2	Immunoprecipitation/Luciferase reporter assay/Western blot	26561204
SGMS1	Immunoprecipitation/Luciferase reporter assay/Western blot	26561204
SAV1	Immunoprecipitation/Luciferase reporter assay/Western blot	26561204
PTEN	Immunofluorescence/Luciferase reporter assay/Western blot	25695913
ST8SIA4	Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	27527856
NDRG2	Luciferase reporter assay/qRT-PCR/Western blot	27533020
TGFBR1	Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	26682928
TGFBR2	Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	26682928
TGFBRAP1	Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	26682928
IL7	ELISA/Luciferase reporter assay/qRT-PCR	25962782
Target genes of hsa-miR-212-3p (Up regulated, log2FC= 2.56, adj P value <0.001)
MECP2	Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot/Reporter assay; Western blot; Microarray;	20020497
TJP1	Western blot/Luciferase reporter assay	18162065
PEA15	qRT-PCR/Luciferase reporter assay/Western blot	20388802
PTCH1	Luciferase reporter assay/qRT-PCR/Western blot	22357618
RB1	Luciferase reporter assay/Western blot/Reporter assay; Western blot; qRT-PCR	21329664
KCNJ2	Immunoblot/Luciferase reporter assay/qRT-PCR	22880819
CCNA2	Western blot; qRT-PCR	21329664
CCNB1	Western blot; qRT-PCR	21329664
RBP2	Luciferase reporter assay/qRT-PCR/Western blot	23922798
MYC	Luciferase reporter assay/qRT-PCR/Western blot	21053104
SOD2	Luciferase reporter assay	23583431
ACHE	In situ hybridization/Luciferase reporter assay/qRT-PCR/Western blot	23974008
IRAK4	Luciferase reporter assay	23264652
PXN	Luciferase reporter assay/qRT-PCR/Western blot	26693054
RFXAP	Luciferase reporter assay/qRT-PCR/Western blot	26337469
SOX4	Immunoprecipitation/Luciferase reporter assay/qRT-PCR/Western blot	26377202
SOX11	Luciferase reporter assay	25766675
ABCG2	Luciferase reporter assay	25766675
ADCY1	Luciferase reporter assay	25766675
SGK3	Luciferase reporter assay/qRT-PCR/Western blot	25720694
SMAD2	GFP reporter assay/qRT-PCR/Western blot	25988335
Target genes of hsa-miR-199a-5p (Up regulated, log2FC= 1.11, adj P value <0.001)
IKBKB	Reporter assay/Western blot/GFP reporter assay/Luciferase reporter assay/qRT-PCR	18408758
LIF	ELISA/Luciferase reporter assay	19011087
JUNB	Luciferase reporter assay/Microarray/qRT-PCR/Western blot	25938468
DDR1	Luciferase reporter assay/qRT-PCR/Western blot	20799954
EDN1	Luciferase reporter assay	19783678
MAP3K11	Luciferase reporter assay/qRT-PCR/Western blot	21048306
HIF1A	Luciferase reporter assay/qRT-PCR/Western blot	21847633
SMARCA2	Luciferase reporter assay/qRT-PCR/Western blot	21189327
CD44	Luciferase reporter assay	22498306
TMEM54	Immunoblot/In situ hybridization/Luciferase reporter assay/Northern blot/qRT-PCR	22903020
SMAD4	Luciferase reporter assay	22821565
SULT1E1	Luciferase reporter assay	23060436
GPR78	Luciferase reporter assay	23085757
ERBB2	Western blot	23437196
UNG	Luciferase reporter assay	23228472
CAV1	Luciferase reporter assay	23459460
SIRT1	qRT-PCR/Western blot	23760629
HSPA5	Luciferase reporter assay/qRT-PCR/Western blot	23598416
ATF6	Luciferase reporter assay/qRT-PCR/Western blot	23598416
ERN1	Luciferase reporter assay/qRT-PCR/Western blot	23598416
KL	EMSA/FACS/Flow/GFP reporter assay/Gluc reporter assay/GUS reporter assay/HITS-CLIP/Immunoblot/Immunocytochemistry/Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot/In situ hybridization	24655788
APOE	ELISA/Luciferase reporter assay/Microarray	23142051
DNAJA4	ELISA/Luciferase reporter assay/Microarray	23142051
ERBB3	Immunoblot/Immunohistochemistry/Luciferase reporter assay/qRT-PCR	23146892
CDH1	Immunofluorescence/Immunohistochemistry/Immunoprecipitation/Luciferase reporter assay/Microarray/qRT-PCR/ Western blot	25080937
PTGS2	Luciferase reporter assay/Western blot	24413338
LIN7A	Luciferase reporter assay	23201090
ARHGAP12	Luciferase reporter assay	23201090
CTSC	Luciferase reporter assay	23201090
RND1	Luciferase reporter assay	23201090
NECTIN1	Luciferase reporter assay	23201090
DRAM1	Luciferase reporter assay/Western blot	23337876
BECN1	Luciferase reporter assay/Western blot	23337876
MAFB	Luciferase reporter assay/Western blot	24391856
WNK1	Luciferase reporter assay/qRT-PCR/Western blot	24299514
NFKB1	Luciferase reporter assay/qRT-PCR/Western blot	24299514
ACVR1B	Immunoprecipitation/qRT-PCR/Western blot/Luciferase reporter assay	25258381
VEGFA	Luciferase reporter assay	24155090
CDH2	Luciferase reporter assay/Western blot	25041364
SNAI1	Luciferase reporter assay/Western blot	25041364
GSK3B	ChIP-seq/FACS/Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot	23764775
FZD4	ChIP-seq/FACS/Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot	23764775
WNT2	ChIP-seq/FACS/Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot	23764775
JAG1	ChIP-seq/FACS/Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot	23764775
HK2	Luciferase reporter assay/qRT-PCR/Western blot	26054020
KRAS	qRT-PCR/Western blot	26027911
SMAD3	Luciferase reporter assay	27038547
ETS1	Luciferase reporter assay	27094578
CCR7	Luciferase reporter assay/qRT-PCR	27814720
PDE4D	Luciferase reporter assay	26344589
TGFB2	Luciferase reporter assay/qRT-PCR	27842518
PIK3CD	Luciferase reporter assay/qRT-PCR	27842518
SETD2	Luciferase reporter assay	26191163
LDLR	Luciferase reporter assay/qRT-PCR/Western blot	26163491
CLTC	Luciferase reporter assay/qRT-PCR/Western blot	26163491
RAB21	Luciferase reporter assay/qRT-PCR/Western blot	26163491
OSCP1	Luciferase reporter assay/Microarray/qRT-PCR/Western blot	26781989
PIAS3	Luciferase reporter assay/qRT-PCR/Western blot	28120918
PSMD9	Luciferase reporter assay/qRT-PCR/Western blot	28120918
CDKN1C	Luciferase reporter assay	25588980
ITGA3	Luciferase reporter assay/qRT-PCR/Western blot	28324890
FZD6	Immunohistochemistry/LacZ reporter assay/qRT-PCR/Western blot	25772759
MAP4K3	Luciferase reporter assay/qRT-PCR/Western blot	28099144
TGFBR1	Luciferase reporter assay/Microarray/qRT-PCR/Western blot	25938468
SLC27A1	Luciferase reporter assay/Microarray/qRT-PCR/Western blot	27562723
Target genes of hsa-miR-339-5p (Up regulated, log2FC= 1.61, adj P value=0.0.3)
BCL6	Immunohistochemistry/Luciferase reporter assay/Western blot	24917186
PTP4A1	Luciferase reporter assay/qRT-PCR/Western blot	23696794
BACE1	Luciferase reporter assay/qRT-PCR/Western blot	24352696
NOVA1	Luciferase reporter assay/qRT-PCR/Western blot	26391641
MDM2	Luciferase reporter assay	24882579
Target genes of hsa-miR-30b-5p (Up regulated, log2FC= −1.12, adj P value<0.001)
CCNE2	Luciferase reporter assay/qRT-PCR	22384020
SMAD1	Luciferase reporter assay/Western blot	22253433
BCL6	Luciferase reporter assay/qRT-PCR/Western blot	22387553
SOCS1	Luciferase reporter assay/qRT-PCR/Western blot	22387553
SNAI1	Luciferase reporter assay/qRT-PCR/Western blot	21633953
CAT	Luciferase reporter assay/qRT-PCR/Western blot	22880027
ATG12	GFP reporter assay/Luciferase reporter assay/qRT-PCR/Western blot	22647547
RUNX2	Microarray/qRT-PCR/Luciferase reporter assay	23316327
DNMT1	qRT-PCR/Western blot	24913034
ERG	Luciferase reporter assay/Microarray/qRT-PCR/Western blot	23728339
BECN1	GFP reporter assay/Luciferase reporter assay/qRT-PCR/Western blot	22647547
SERPINE1	GFP reporter assay/Luciferase reporter assay/Northern blot/qRT-PCR/Western blot	25170877
DLL4	In situ hybridization/Luciferase reporter assay/qRT-PCR/Western blot	23086751
SIX1	Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	24593661
BCL9	Immunofluorescence/Immunohistochemistry/In situ hybridization/Luciferase reporter assay/qRT-PCR/Western blot	24599134
NOTCH1	Luciferase reporter assay/qRT-PCR/Western blot	24029422
MBNL1	Luciferase reporter assay	26840300
MBNL2	Luciferase reporter assay	26840300
PDGFRB	Luciferase reporter assay/qRT-PCR/Western blot	23893664
TP53	Luciferase reporter assay/qRT-PCR/Western blot	24029422
BCL2	Luciferase reporter assay	24178239
MBNL3	Luciferase reporter assay	26840300
EIF2S1	Luciferase reporter assay	26898246
EIF5A2	Luciferase reporter assay/qRT-PCR/Western blot	26309359
RASAL2	Luciferase reporter assay	25547110
CTHRC1	Luciferase reporter assay/Western blot/qRT-PCR	26388700
HOXA1	Luciferase reporter assay/qRT-PCR/Western blot	28189678
Target genes of hsa-miR-149-5p (Down regulated, log2FC= −1.44, adj P value <0.001)
SP1	Luciferase reporter assay	22821729
FOXM1	Luciferase reporter assay/Western blot	23762558
ZBTB2	Immunohistochemistry/Luciferase reporter assay/qRT-PCR/Western blot	23144691
GIT1	Immunofluorescence/Immunohistochemistry/Luciferase reporter assay/Microarray/qRT-PCR/Western blot	24608434
FGFR1	Luciferase reporter assay/qRT-PCR/Western blot	24463821
GPC1	Luciferase reporter assay/qRT-PCR/Western blot	24463821
MYD88	ELISA/qRT-PCR/Western blot	24375488
IL6	ELISA/Luciferase reporter assay/qRT-PCR/Western blot	25916550
PTGER2	ELISA/Luciferase reporter assay/qRT-PCR/Western blot	25916550
BBC3	Immunofluorescence/Luciferase reporter assay/Northern blot/qRT-PCR/Western blot	23873935
FGF21	Western blot	27061435
MTHFR	Immunoblot/Luciferase reporter assay/qRT-PCR/Western blot	28045918
FASLG	Luciferase reporter assay/qRT-PCR/Western blot	28013316
PPM1F	Luciferase reporter assay/Western blot	26498692

Protein Interaction Network Construction

In order to further understand the biological interaction of SIX1, CCNA2, CCNB1, FOXM1, and RUNX2, we constructed an expanded PPI network. The network had 55 nodes, 52 genes, and 951 edges, with an average local clustering coefficient of 0.836 (Figure 3). These 52 genes represent a gene cluster based on SIX1, CCNA2, CCNB1, FOXM1, and RUNX2.

Figure 3

A PPI network with 55 nodes, 52 genes, and 951 edges was constructed based on five DE-miRNAs’ target mRNAs (SIX1, CCNA2, CCNB1, FOXM1, and RUNX2).

Expression Analyses of Target mRNAs of DE-miRNAs

Next, in order to investigate the biological role of this gene cluster in keloids, we performed gene enrichment analysis. DAVID online tool was used for GO annotation and KEGG pathway analyses. The top five items in each enrichment analysis are shown in Figure 4A. For the purpose of validating the expression of those five genes in keloid, we further analyzed the expression of target genes in the single-cell type clusters identified in skin, and the data was visualized by a bar chart (Figure 4B). The single-cell RNA sequencing (scRNAseq) data were obtained from the HPA database. The results showed that the expression level of CCNA2, CCNB1, RUNX2, and FOXM1 was relatively higher in fibroblast c-9 cell cluster (fibroblast mesenchymal cells) than in other cell clusters in the skin. SIX1 expression increased in the fibroblast c-3 cell cluster, which also identified as fibroblast mesenchymal cells. The expression profile of the five genes in a specific tissue or organ system was also identified using BioGPS. The results showed that CCNA2 expression was highest in B lymphoblast, followed by CD34 and CD105+ endothelial cells. Similarly, CCNB1 expression was highest in B lymphoblast, followed by CD71+ early erythroid and CD105+ endothelial cells. FOXM1 expression was highest in Burkitt’s lymphoma, followed by CD71+ early erythroid and B lymphoblast. SIX1 and RUNX2 were not specifically expressed in tissues. Next, we tried to analyzed the expression and existence of these five genes in healthy human tissues and organs through Bgee database. The data were expressed in a bar chart (Figure 4C). We found that compared with normal subcutaneous adipose tissue, except SIX1, the expression difference of other genes in subcutaneous adipose tissue and skin was not obvious. This suggests that the expression differences of CCNA2, CCNB1, RUNX2, and FOXM1 in tissue adjacent to the skin are not obvious in healthy donors. Therefore, we suspect that abnormal expression of CCNA2, CCNB1, RUNX2, and FOXM1 may be associated with keloid formation.

Figure 4

(A) The enrichment analysis of the 52 genes mentioned in Figure 3 was analyzed to study the biological function of the gene cluster based on five DE-miRNAs’ target mRNAs. (B) We further analyzed the expression of five target genes in the single-cell-type clusters identified in the skin to analyze the expression of those five mRNAs in different cell types. (C) Expression and existence of these five genes in different healthy human skin tissues was analyzed using Bgee database.

Validation of the Identified miRNAs/mRNAs

After we screened a gene cluster containing five mRNAs and three miRNAs through the above steps, for further confirmation, the expression level of five genes (SIX1, CCNA2, CCNB1, FOXM1, RUNX2) and three miRNAs (miR-30b-5p, miR-212-3p, miR-149-5p) in keloid tissue and normal skin was analyzed by using RT-PCR. According to our results, we found that the expression of miR-30b-5p, miR-212-3p, CCNB1, RUNX2 and SIX1 in keloid was significantly higher than that in healthy donors, and the expression of miR-149-5p in keloid was significantly lower than that in healthy donors, and the difference was statistically significant (Figure 5).

Figure 5

The expression of DE-miRNAs and DE-target mRNAs was validated using RT-PCR. The 2-ΔΔCt method was used to determine the gene expression levels. Differences between groups were analyzed using Student’s t-test. The significance was defined at P < 0.001 (***), P < 0.01 (**), and P < 0.05 (*).

Discussion

Herein, this study aimed to identify biomarkers for keloid and a cluster containing five genes (SIX1, CCNA2 CCNB1, FOXM1, RUNX2) controlled by three miRNAs (miR-30b-5p, miR-212-3p, miR-149-5p) were determined. In this study, we identified a novel pathway, miR-149-5p-FOXM1. miR-149-5p is overexpressed in multiple cancers, such as breast cancer, lung cancer, kidney cancer, and bladder cancer. By regulating il-6/STAT3 pathway, miR-149-5p thus affects the proliferation and metabolism of tumor cells.13–17 However, there has been no research report on the correlation between the expression of miR-149-5p and keloid so far. This is the first study to reveal the association between miR-149-5p and keloid. FOXM1 (Forkhead Box M1) is involved in cell proliferation. FOXM1 protein is phosphorylated in the M phase of the cell cycle, thus regulating the expression of several cell cycle genes, such as cyclin B1 and cyclin D1. Consistent with previous studies, enrichment analysis in this study demonstrated that the PPI network containing FOXM1 is mainly associated with the cell cycle. Previous research showed that FOXM1 knockout inhibits keloid fibroblast activation and extracellular matrix production by affecting the activation of the TGF-β1/Smad pathway.18 Another research also found that FOXM1 promotes keloid progression by participating in the HOXA11-AS-miR-205-5p-FOXM1 pathway.19 Therefore, the novel pathway miR-149-5p-FOXM1 could provide the underlying mechanism of keloid.

Our results indicate that miR-212-3p can promote keloid formation by regulating CCNA2 or CCNB1. miR-212-3p regulates cancer tumorigenesis, differentiation, proliferation, and survival by inhibiting major genes, including SGK3, SOX4, SOX11, SOX5.20–24 There are currently no reports in literature for the function of miR-212-3p on keloid. Unlike miR-212-3p, CCNB1 and CCNB2 alterations may contribute to keloid formation via dysregulation of cell cycle and p53 signaling pathways.25 Furthermore, another research showed that CCNB1 expression might be related to 5-fluorouracil treatment in keloid fibroblasts, suggesting that further study of CCNB1 may lead to new breakthroughs in the treatment of keloids.26 In our study, we not only confirmed the abnormal expression of miR-212-3p and CCNB1 in keloid compared with normal skin but also found that the expression differences of CCNB1 and CCNA2 in normal tissue adjacent to the skin are not obvious in healthy people. Therefore, we suspect that abnormal expression of miR-212-3p and CCNB1 may play a role in keloid formation. Considering strong biological connections between CCNB1, CCNB2, and CCNA2, they may act together to influence the keloid formation via regulating fibroblast cell cycle, and miR-212-3p may act as a key regulatory node in this process.

Our findings suggest that the miR-30b-5p-SIX1/RUNX2 pathway may be associated with keloid. Various studies have shown that miR-30b-5p is associated with tumorigenesis and tumor metastasis.27 Previous studies have illustrated that miR-30b-5p promotes keloid fibroblast proliferation by participating in DLEU2-miR-30a-5p/miR-30b-5p pathway.28 RUNX2 is a member of the RUNX family of transcription factors and encodes a nuclear protein with a Runt DNA-binding domain. Consistent with our study, previous studies have indicated that keloid is related to the expression of RUNX2.29 Specifically, RUNX2 siRNA therapy may reverse pathological keloid phenotypes by inhibiting PI3K/AKT signaling,30 suggesting that RUNX2 is a promising biomarker for keloid therapy. Besides, previous studies have shown that miR-194-3p can directly regulate RUNX2.31 Herein, RUNX2 was associated with miR-30b-5p, providing a new regulatory target for RUNX2 regulation in keloid. The SIXI gene located on chromosome 14 has been associated with limb development, and it encodes a homeobox protein called SIX protein. The autosomal dominant deafness type 23 and branchiootic syndrome type 3 are associated with SIX1. However, no study has comprehensively demonstrated the correlation between SIX1 and keloid. In this study, SIX1 was differentially expressed in KF vs CF, suggesting that SIX1 may promote the formation of the keratinocyte layer of keloid. The abnormal expression of miR-30b-5p, RUNX2 and SIX1 in keloid compared with normal skin was validated. In addition, we found that there was no obvious difference in RUNX2 expression between normal skin tissue and adjacent adipose tissue, suggesting that the differential expression of RUNX2 and SIX1 in keloid may be related to the pathogenesis of keloid, which may be regulated by Mir-30b-5p.

However, this study has several limitations. First, the small number of samples obtained from GSE113620 introduced bias during the analysis of DE-miRNAs. Therefore, more chip data are necessary for statistical analysis to verify the experimental results. Moreover, the functional and molecular mechanisms of genes are complex, and predictions based on bioinformatics alone require cell and animal experiments. Lastly, in the three RNA-seq analysis (GSE113620, GSE44270, GSE92566) we studied, the source of sample was different. In GSE44270, sample were excised from patients undergoing surgery, and then primary cells of skin and scar tissue were cultured for further microarray analysis. However, in GSE113620 and GSE92566, skin and scar tissues were excised and preserved with RNAlater, and subsequently the sample tissue was used for RNA extraction and RNA-seq analysis. We believe that the formation of keloid is closely related to its environment, such as tension and inflammation. Therefore, after the removal of external factors, the gene expression in cultured primary cell may be different from that of tissue sample. That is to say, the results of this study need to be further verified with appropriate sample.

Conclusion

In our study, a PPI network containing 52 genes was constructed based on five target mRNAs of DE-miRNAs. Gene enrichment analysis showed that the biological function of this gene cluster is mainly associated with the cell cycle. Therefore, the three DE-miRNAs may influence keloid proliferation by affecting the gene functional groups centered on CCNA2, CCNB1, SIX1, RUNX2, and FOXM1. This study may provide reliable biomarkers for accurate diagnosis and individualized keloid treatment, and further clinical application.