Identifying miRNAs Associated with the Progression of Keloid through mRNA-miRNA Network Analysis and Validating the Targets of miR-29a-3p in Keloid Fibroblasts

Background Keloid has brought great trouble to people and currently has no uniformly successful treatment. It is urgent to find new targets to effectively prevent the progress of keloid. The current research mainly identifies the differentially expressed genes (DEGs) in keloid through high-throughput sequencing technology and bioinformatics analysis technology, to screen new therapeutic targets and potential biomarkers. However, due to the different samples, different control groups, and small sample sizes, the sequencing results obtained from different studies are quite different and lack reliability. It is necessary to analyze the existing datasets in a reasonable way. Methods Datasets about keloid were filtered in Gene Expression Omnibus (GEO) and ArrayExpress databases according to the inclusion and exclusion criteria. The discovery datasets were used for summarizing significant DEGs, and the validation datasets were to validate the mRNA and miRNA expression levels. The Encyclopedia of RNA Interactomes (ENCORI) online platform was used to predict the interactions between miRNAs and their target mRNAs. Protein-protein interaction network (PPI network) analysis and functional enrichment analysis were conducted. miRNA-mRNA network was established by Cytoscape software and verified in keloid tissue (n = 8) by RT-qPCR. miR-29a-3p mimic and inhibitor were transfected into keloid fibroblasts (KFs) to preliminary verify its targets, the prognostic value of which was estimated by the receiver operating characteristic (ROC) curve. Results A total of 6 datasets involving 20 patients were included. 15 miRNAs and 12 target mRNAs were identified as potential biomarkers for keloid patients. The RT-qPCR results showed that miR-29a-3p, miR-92a-3p, and miR-143-3p were downregulated, and all their target mRNAs were upregulated in keloid tissue (P < 0.05). The expression of COL1A1, COL1A2, COL3A1, COL5A1, and COL5A2 decreased when miR-29a-3p was overexpressed but increased when miR-29a-3p was knocked down (P < 0.05). And these genes had a good performance in the diagnosis of keloid, especially when using keloid nonlesional skin or normal scar tissues as controls. Conclusion The miRNA-mRNA network, especially miR-29a-3p and its targets, may provide insights into the underlying pathogenesis of keloid and serve as potential biomarkers for keloid treatment.


Introduction
Keloid is a complex fibroproliferative disease similar to the benign tumor growing outside the edge of the original wound, which can cause skin dysfunction and aesthetic deformity by invading adjacent normal tissues [1]. In addition, a cross-sectional survey on the burden of keloid shows that keloid has caused serious damage to the emotional health of patients [2]. In a word, keloids pose a great threat to people's physical and mental health, but the underlying pathogenesis of keloid formation remains to be elucidated, and the treatments for keloid are not satisfactory (including triamcinolone acetonide injection, surgical resection, and radiotherapy) [3,4]. Therefore, it is urgent to seek effective biomarkers of keloid, which can further reveal the molecular mechanism of keloid formation and provide new insights for the development of treatment methods.
miRNA is a kind of noncoding single-stranded RNA with a length of about 22 nucleotides. It can bind hundreds of target genes and regulate their expression. A large number of studies have confirmed that the miRNA-mRNA axis affects the development of keloid by regulating cell proliferation, apoptosis, migration, invasion, and collagen synthesis [5]. Therefore, miRNAs that play an important role in the development of keloid may become new targets for the treatment of keloid [6].
The current research identifies abnormally expressed genes to screen new therapeutic targets and potential biomarkers mainly through high-throughput sequencing technology and bioinformatics analysis technology [7][8][9]. However, the sequencing results obtained from different studies are quite different. There are three main reasons: (1) the high cost of high-throughput sequencing technology makes it difficult to sequence large-scale samples. (2) The samples used for detecting are different (including keloid tissue, keloid fibroblast, and keloid keratinocyte).
(3) The source of normal skin used as control is different (including normal skin of keloid patients, normal skin of healthy people, and foreskin). Therefore, to obtain reliable results, it is necessary to analyze the existing datasets in a reasonable way.
In our study, datasets about keloid were filtered in Gene Expression Omnibus (GEO) and ArrayExpress databases according to the inclusion and exclusion criteria. Four datasets (GSE92566, GSE83286, GSE158395, and GSE190626) were screened and treated as discovery datasets for summarizing significantly differentially expressed genes (DEGs), and the other two datasets (GSE113619 and GSE113620) were served as validation datasets to validate the mRNA and miRNA expression levels. In addition, the functions of these DEGs were annotated by Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Finally, a miRNA-mRNA network was established and verified in keloid tissue, and the targets of miR-29a-3p were preliminary verified in keloid fibroblasts (KFs). Our research was aimed at finding novel therapeutic targets for keloid.

Clinical Specimen Information and Fibroblast Culture.
Lesional and nonlesional keloid skin samples were collected from patients undergoing keloid surgery in the Xiang'an Hospital affiliated with Xiamen University. The clinical specimen information is presented in Table 1. This experiment was conducted in accordance with the Declaration of Helsinki and with the approval of the Medical Ethics Committee of Xiang'an Hospital Affiliated with Xiamen University and informed patient consent in writing.
The primary fibroblasts were successfully isolated based on the previously reported procedures. Briefly, we removed the epidermis and fat of the skin specimen. The dermal tissue was cut into 5 mm 2 explants and was evenly distributed in the sterile Petri dish. The explants were incubated in Dulbecco's Modified Eagle Medium supplemented with 10% fetal bovine serum and 2% penicillin-streptomycin in an atmosphere of 5% carbon dioxide for the primary and subsequent cultures. The KFs at passages 3-5 were carried out to subsequent experiments.

2.2.
Datasets. The Gene Expression Omnibus (GEO) (https:// www.ncbi.nlm.nih.gov/geo/) and ArrayExpress (https://www .ebi.ac.uk/arrayexpress/) databases were searched to identify relevant transcriptomic profiling datasets. Datasets containing keloid transcriptomic profiling in Homo sapiens were potentially eligible. Datasets were excluded if they were not assaying keloid tissue samples, not measuring mRNA or miRNA, and not taking the nonlesional skin of keloid patients themselves as the control group. Furthermore, datasets that did not contain complete expression data were also excluded. Additional datasets could be added by manual search of the reference of included studies. The data of this part is obtained from the public database, so it is not necessary to obtain the patient's informed consent and ethical review.
2.3. Data Processing. R software (version 4.0.3) was used to convert the ID of datasets, and probes with missing gene symbols were excluded. The "LIMMA" software package was used to analyze the DEGs in the datasets (adj:P:value < 0:05, jlogFCj > 1). These DEGs were visualized by volcano plots. We screened common differentially expressed genes (co-DEGs) existing in at least 3 discovery datasets and divided co-DEGs into two groups, namely, co-up DEGs and codown DEGs.

Functional Enrichment
Analysis. Gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of co-up DEGs and co-down DEGs are carried out by using Bioinformatics online platform (http://www.bioinformatics.com.cn). GO enrichment analysis includes biological processes (BP), molecular functions (MF), and cellular components (CC).

Protein-Protein Interaction (PPI) Network Analysis and
Identification of Hub Genes. STRING online platform (https://cn.string-db.org/) was used to construct the PPI network which was analyzed by Cytohubba (a plug-in of Cytoscape). The top 20 genes ranked by different methods were  Table 2. U6 and GAPDH were used as internal controls. The fold change in the relative gene expression to the con-trol levels was determined using the standard 2−ΔΔCt method.

Cell Transfection.
Mimics or inhibitors targeting miR-29a-3p and their respective blank controls were designed by Ruibo (Wuhan, China). To achieve better transfection efficiency, primary fibroblasts were seeded on 6-well plates in a conditioned medium without antibiotics and serum at a density of 5 × 10 4 cells/mL. INTERFERin® (Polyplus, France) was applied to transfect these miRNAs into KFs according to the manufacturer's protocol. After transfection for 48 h, the expression of miR-29a-3p and the target genes were detected through qPCR.
2.9. Evaluation of Diagnostic Value. Different datasets were used to evaluate the diagnostic value of key genes. GSE92566, GSE83286, GSE158395, and GSE190626 were used to evaluate the diagnostic value of key genes between keloid lesional skin and nonlesional skin; GSE113619 was used to evaluate the diagnostic value of key genes between the normal skin from keloid-prone or normal individuals; GSE188952 was used to evaluate the diagnostic value of key genes between keloid tissue and hypertrophic scar or normal scar.

Statistical Analysis.
All experiments were performed at least three times. The data are presented as the mean ± Reverse 5′-GCTGGTGGTCCAGGGGTCTTACT-3′ 3 BioMed Research International standard deviation (S.D.). Nominal variables were expressed as a count or as a percentage. The statistical comparisons between the groups were performed by paired Student's t -test or unpaired Student's t-test or one-way analysis of variance (ANOVA) using GraphPad Prism V. 9 (GraphPad Software, San Diego, CA, USA). P < 0:05 indicated a significant difference.

Datasets Search and Characteristics of Included Datasets.
A detailed overview of the selection process is shown in Figure 1. According to the inclusive and exclusive criteria, a total of 6 studies, including 20 patients, were enrolled in the analysis. In addition, a dataset (GSE188952) using hypertrophic scars and normal scars tissue as control was also included in the study to further evaluate the diagnostic value of key genes in keloid. The characteristics of the datasets included in the analysis are described in Table 3.  Table 4.

Functional Enrichment
Analysis. 318 co-up DEGs and 175 co-down DEGs were analyzed by GO and KEGG analysis to explore their biological functions. It is worth noting that the 318 co-up DEGs in keloid are mainly related to the organization of extracellular matrix, extracellular structure, and collagen fibril, also involved in the development of bone, cartilage, and connective tissue (Figure 4(a)). The 175 co-down DEGs are mostly involved in sodium and potassium ion homeostasis and skin development ( Figure 4(b)). In addition, the co-up DEGs and co-down DEGs are closely related to the focal adhesion ( Figure 5(a)) and peroxisome proliferator-activated receptor (PPAR) signal pathways ( Figure 5(b)) respectfully. A previous study had proved that increased transcriptional response to mechanical strain in keloid fibroblasts due to increased focal adhesion complex formation [11]. And Zhang et al. reported that troglitazone (PPAR-γ agonist) depressed transforming growth factor-β1-(TGF-β1-) stimulated collagen type I expression and collagen synthesis in keloid fibroblasts in a concentration-dependent manner [12]. These studies further confirmed the association between the focal adhesion and PPAR signal pathway and keloid.

BioMed Research International
The expression level of the 16 genes was verified by the GSE113619 dataset. The results showed that the expression of 12 genes (Hub genes) significantly increased in the development of keloid (P < 0:05) while that of TGFB3, SOX9, and BGN did not change significantly (P > 0:05) (Figure 7). Unexpectedly, in the discovery datasets, the expression of FN1 was significantly upregulated in keloid but significantly downregulated in the GSE113619 dataset ( Figure 7). Considering that the samples of GSE113619 were collected from the early stage of keloid formation and the samples of discovery datasets were mature keloid, this difference may be caused by different stages of keloid development, which needs more in-depth study to explore the characteristics at different stages of keloid.       (Figure 8). We finally formed a miRNA-mRNA network composed of 12 mRNAs and 15 miRNAs by Cytoscape (Figure 9).  To further verify the miRNA-mRNA network, we detected the lesional and nonlesional keloid skin samples of 8 keloid patients by RT-qPCR. Considering that all Hub genes were upregulated and miRNA usually inhibits the expression of target genes, we mainly verified the 3 downregulated miRNAs and their upregulated targets. The results of RT-qPCR showed that miR-29a-3p, miR-92a-3p, and miR-143-3p showed a significant downregulation trend in keloid tissue, and the increase of their targets was also significant ( Figure 10).

miR-29a-3p
Regulated the Expression of COL1A1, COL1A2, COL3A1, COL5A1, and COL5A2 in KFs. The above results showed that miR-29a-3p was significantly downregulated in the development of keloid. More importantly, the miRNA-mRNA network showed that it can regulate the expression of 9 Hub genes (upregulated). Therefore, we speculated that miR-29a-3p plays a central role in Hub miRNA. Previous studies have confirmed that miR-29a-3p was significantly downregulated in keloid tissues and fibroblasts and can regulate the proliferation, apoptosis, migration, and invasion of KFs by targeting COL1A1 and COL3A1 [13,14]. However, more potential targets for the role of miR-29a-3p remain to be clarified. To further verify the regulation of miR-29a-3p on the expression of these 9 Hub genes, we transfected miR-29a-3p mimic and inhibitor into KFs to detect the expression level of these Hub genes. RT-qPCR results showed that the expression of COL1A1, COL1A2, COL3A1, COL5A1, and COL5A2 decreased when miR-29a-3p was overexpressed but increased when miR-29a-3p was knocked down, indicating that these Hub genes may be potential targets for miR-29a-3p to regulate the biological behavior of KFs ( Figure 11).

Discussion
Keloid has brought great trouble to people and currently has no uniformly successful treatment [3]. Therefore, it is urgent to find new targets to effectively prevent the progress of keloid. In this study, a miRNA-mRNA network closely related to keloid was constructed through a comprehensive analysis of multiple datasets. miR-29a-3p was considered to be the real central miRNA through bioinformatics and molecular biology experiments. It was downregulated in keloid and can negatively regulate the expression of COL1A1, COL1A2, COL3A1, COL5A1, and COL5A2. It is suggested that miR-29a-3p plays an important role in the development of keloid and may become an effective target for the treatment of keloid.
The dysregulated expression of miR-29a-3p plays a role in a variety of fibrotic diseases including keloid. Zhang et al. found that the expression of miR-29a was significantly decreased in keloid tissues and fibroblasts (n = 9) compared with healthy controls (n = 6). They also found that the mRNA and protein levels of type I and type III collagen decreased in KFs with miR-29a overexpression and confirmed that it can directly regulate the expression of COL3A1 by luciferase gene reporting experiment [13]. Wang et al. confirmed that the expression of miR-29a in keloid tissue and fibroblasts (n = 80) was significantly lower than that in normal skin (n = 91) and could inhibit the viability, migration, invasion, and EMT of KFs while promoting its apoptosis. In addition, they proposed that miR-29a plays the above functions mainly by directly inhibiting the expression of COL1A1 [14]. Wu et al. found that miR-29a-3p was upregulated in hypertrophic scar (HS) tissue and fibroblasts compared with normal skin tissue and confirmed that miR-29a-3p inhibited the proliferation of hypertrophic scar fibroblasts (HSF) and the production of extracellular matrix (ECM) by directly targeting FRS2 [15]. Yuan et al. found that the expression level of miR-29a in mouse scar tissue was significantly downregulated compared with normal skin tissue and proposed that exosomes from miR-29a modified adipose mesenchymal stem cells could reduce excessive scar formation by inhibiting TGF-β2/Smad3 signal [16]. Nevertheless, more potential targets of miR-29a-3p have yet to be clarified.

12
BioMed Research International type I collagen and COL3A1 encodes type III collagen, the major components of the extracellular matrix in keloid tissue. Compared with normal skin, type I and III collagen deposition in keloid tissue increased significantly, and previ-ous studies have confirmed that miR-29a-3p can regulate the biological behavior of KFs by directly targeting COL1A1 [14] and COL3A1 [13]. COL5A1 and COL5A2 are involved in encoding type V collagen, and its mutation is related to

BioMed Research International
Ehlers-Danlos syndrome. Yokota et al. found that the expression of COL5A1 and COL5A2 was significantly increased in the early stage after acute ischemic heart injury, while the decrease of type V collagen leads to the contradictory increase of scar size after infarction [17]. However, there is no research report on the correlation between COL5A1 and COL5A2 expression and keloid so far.
In addition, Corrie et al. conducted a double-blinded, placebo-randomized, within-subject controlled clinical trial and found that remlarsen (a miR-29 mimic) repressed collagen expression and the development of fibroplasia in inci-sional skin wounds. They suggested that remlarsen may be an effective therapeutic to prevent the formation of a fibrotic scar (hypertrophic scar or keloid) or to prevent cutaneous fibrosis (such as scleroderma), which also provides strong human experimental evidence for the clinical application of miR-29a-3p in the treatment of keloid [18].
Except for miR-29a-3p, we also verified the expression level of 2 downregulated miRNA (miR-92a-3p and miR-143-3p). Mu et al. found that miR-143-3p was significantly downregulated in HS tissues and fibroblasts and confirmed that miR-143-3p inhibits HS formation by regulating the  proliferation and apoptosis of human HSF [19]. Wei et al. found that miR-143-3p was downregulated in HS compared with normal skin tissue, and the overexpression of miR-143-3p inhibited the proliferation and invasion of HSF and promoted its apoptosis. They also confirmed that miR-143-3p played the above role by directly targeting Smad3 [20].
Besides, miR-92a-3p [21] was also involved in the development of various fibrotic diseases, but the correlation with keloid has not been reported.
Other target genes in the miRNA-mRNA network were also potentially involved in the development of keloid. For example, COL6A2 and COL6A3 jointly encode type VI  Figure 12: Diagnostic value of COL1A1, COL1A2, COL3A1, COL5A1, and COL5A2 in keloid. 16 BioMed Research International collagen, and their mutation will lead to an autosomal dominant disease, namely Bethlem myopathy. James, Aralikatte, Edmar, Constanza et al. reported the formation of spontaneous keloid in patients with Bethlem myopathy, suggesting that there is a potential correlation between the imbalance of type VI collagen and the occurrence of keloid [22][23][24][25].
In addition, Georgios et al. found that type VI collagen was expressed in the extended edge of human keloid samples and confirmed that type VI collagen is a key regulator of dermal matrix assembly, composition, and fibroblast behavior and may play an important role in wound healing and tissue regeneration [26]. Periostin (POSTN) is a matricellular protein normally expressed in adult skin. Previous studies reported that POSTN was overexpressed in keloid and hypertrophic scars and increased POSTIN expression affects the proliferation, collagen synthesis, migration, and invasion of keloid fibroblasts under hypoxic conditions [27,28]. In a word, the miRNA-mRNA network we constructed intuitively presents the regulatory effect of Hub miRNAs on Hub genes, effectively reveal the key genes and regulatory targets in the development of keloid, and lay the foundation for the development of new treatment schemes.
The current research mainly identifies the DEGs in keloid through high-throughput sequencing technology and bioinformatics analysis technology, to screen new therapeutic targets and potential biomarkers [7][8][9]. However, due to the different samples, different control groups, and small sample sizes, the sequencing results obtained from different studies are quite different and lack reliability. In this study, 6 datasets consisting of 20 keloid patients were screened from the GEO and ArrayExpress databases according to the inclusion and exclusion criteria. The samples and control groups were strictly controlled, and the sample size was expanded to a certain extent. GSE113619 and GSE113620 datasets were used as validation datasets to verify the expression levels of mRNA and miRNA, respectively, and RT-qPCR was used to verify the miRNA-mRNA network, which enhanced the reliability of the analysis results and reduced the false positive rate.
However, our research still has the following limitations: (1) there was no classified analysis of keloid patients with different races, different causes, and different parts; (2) focus on obtaining the most important genes and lack of attention to other genes; (3) the regulation of target genes by miRNA was predicted by ENCORI and has not been further verified by experiments; and (4) the potential target of miR-29a-3p has only been preliminarily verified, and more in-depth research is needed in in vivo and in vitro models.

Conclusions
In conclusion, DEGs between keloid and normal skin tissue were identified and functionally annotated through the analysis of gene expression data from multiple datasets. The miRNA-mRNA network was constructed, and the real central gene miR-29a-3p was screened. The potential target of miR-29a-3p was preliminarily verified, which provides a potential treatment option for the treatment of keloid.