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

Introduction

Keloid is excessive fibrosis of the skin that extends beyond the area of injury and does not regress.¹ Keloid can occur in the joints and mouth after several years of severe injury, including burns, chemical injury, wound, and surgical incision.² Keloids on the joints affect the quality of life, physical, mental, and social health.³ 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.⁴ 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.⁷ Several studies have shown that miRNAs are biomarkers for keloid.⁸ 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; http://www.ncbi.nlm.nih.gov/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 (http://www.bioconductor.org/packages/2.8/bioc/html/affy.html) (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.⁹ The common mRNAs between the target mRNAs and the DEGs in GSE44270 and GSE92566 were then identified. An online Venn tool (http://jvenn.toulouse.inra.fr/app/index.html) was used for Venn diagram analysis.¹⁰

Protein–Protein Interaction (PPI) Network Analysis

A String database (https://string-db.org/) (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 (http://biogps.org, 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 https://bgee.org/ and the development pack BgeeDB.

Cell-Specific Gene Expression Analysis

HPA (http://www.proteinatlas.org/) 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.¹¹ 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.¹¹ 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.¹¹

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.¹² 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

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

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.

Table 3 Target Genes of Identified DE-miRNAs

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.

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.¹⁸ Another research also found that FOXM1 promotes keloid progression by participating in the HOXA11-AS-miR-205-5p-FOXM1 pathway.¹⁹ 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.²⁵ 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.²⁶ 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.²⁷ Previous studies have illustrated that miR-30b-5p promotes keloid fibroblast proliferation by participating in DLEU2-miR-30a-5p/miR-30b-5p pathway.²⁸ 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.²⁹ Specifically, RUNX2 siRNA therapy may reverse pathological keloid phenotypes by inhibiting PI3K/AKT signaling,³⁰ suggesting that RUNX2 is a promising biomarker for keloid therapy. Besides, previous studies have shown that miR-194-3p can directly regulate RUNX2.³¹ 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.