Identification of crucial noncoding RNAs and mRNAs in hypertrophic scars via RNA sequencing

Hypertrophic scarring (HS) is a dermal fibroproliferative disorder characterized by excessive deposition of collagen and other extracellular matrix components. The aim of this study is to explore crucial long noncoding RNAs (lncRNAs) and circular RNAs (circRNAs) associated with HS and provide a better understanding of the molecular mechanism of HS. To investigate the lncRNA, circRNA and mRNA expression profiles, we performed RNA sequencing of human HS and normal skin tissues. After the identification of differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs) and circRNAs (DEcircRNAs), we performed functional enrichment of DEmRNAs. Further on, we constructed DElncRNA/DEcircRNA–DEmRNA coexpression networks and competing endogenous RNA regulatory networks, and performed functional analyses of the DEmRNAs in the constructed networks. In total, 487 DEmRNAs, 92 DElncRNAs and 17 DEcircRNAs were identified. DEmRNAs were significantly enriched in processes such as collagen fibril organization, extracellular matrix–receptor interaction and the phosphatidylinositol 3‐kinase (PI3K)–Akt signaling pathway. In addition, we detected 580 DElncRNA–DEmRNA and 505 DEcircRNA–DEmRNA coexpression pairs. The competing endogenous RNA network contained 18 circRNA–microRNA (miRNA) pairs, 18 lncRNA–miRNA pairs and 409 miRNA–mRNA pairs, including 10 circRNAs, 5 lncRNAs, 15 miRNAs and 160 mRNAs. We concluded that MIR503HG/hsa‐miR‐204‐3p/ACAN, MIR503HG/hsa‐miR‐431‐5p/TNFRSF9, MEG3/hsa‐miR‐6884‐5p/ADAMTS14, AC000035.1‐ADAMTS14 and hsa_circ_0069865‐COMP/ADAM12 interaction pairs may play a central role in HS.

Scar formation is an inevitable result of wound healing. Hypertrophic scarring (HS), a type of pathological scarring, with a protruding surface, irregular shape and burning and itching sensations on the skin surface, often significantly affects patients' quality of life [1]. HS is characterized by excessive deposition and alterations in morphology of collagen and other extracellular matrix (ECM) proteins [2]. Clinically, it is identified by excessive dermal fibrosis and scarring resulting from the imbalance between collagen synthesis and degradation during wound healing [3]. Although numerous interventions for HS, including surgical removal, radiotherapy, steroid injection and cryotherapy, are available, these treatments cannot achieve a stable curative effect [4]. The etiology and pathogenesis of HS have been explored for decades, but the molecular mechanisms of HS remain poorly understood [5]. Therefore, it is of great importance to elucidate the mechanisms of HS and explore new targets for treatment of HS.
Long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), as novel noncoding RNAs, have been reported to be involved in HS. For instance, Li et al. [6] indicated that up-regulated lncRNA8975-1 in HS fibroblasts inhibited fibroblast proliferation and reduced collagen expression. Nong et al. [7] demonstrated that lncRNA COL1A2-AS1 inhibited fibroblasts proliferation to suppress HS formation via regulating the miR-21-Smad7 pathway. Wu et al. [8] reported that miR-155 inhibited the formation of HS fibroblasts by targeting hypoxia induciblefactor 1 subunit alpha (HIF-1a) via the PI3K-AKT pathway. Shen et al. [9] suggested that miR-145-5p arrested the development of fibrogenesis and decreased HS formation by reducing the expression of Smad2/3. Zhang et al. [10] found that miR-137 inhibited proliferation and metastasis of HS fibroblasts via targeting pleiotrophin. However, to the best of our knowledge, there was only one study exploring the expression profiles of circular RNA (circRNA) in HS [11].
This study investigated the mRNA, lncRNA and circRNA expression profiles of HS to identify the differentially expressed mRNAs (DEmRNAs), lncRNAs (DElncRNAs) and circRNAs (DEcircRNAs) associated with HS. In addition, a competing endogenous RNA (ceRNA) (DEcircRNA/lncRNA-miRNA-DEmRNA) regulatory network was conducted. This study sought to make a contribution to elucidating the potential molecular mechanisms of HS and lay a foundation for the treatment of HS.

Subjects and samples
HSs and adjacent normal skin tissues were collected from three male patients in our hospital, aged 6, 21 and 22 years. The site of HS was neck, left arm and right arm. All samples were collected after obtaining written informed consent from every participant. This study was approved by the ethics committee of The 980st Hospital of the PLA Joint Logistics Support Force (2020-KY-25) and performed in accordance with the Declaration of Helsinki. Total RNAs were isolated from HS and normal skin tissues with TRIzol reagent. Based on the Illumina HiSeq X-ten platform, sequencing was performed.

Quality control of raw sequencing and mapping of clean reads
To obtain clean reads from RNA sequencing results, we removed sequences with low quality, including adapter sequences, sequences with quality score < 20, sequences with N base rate of raw reads > 10% and sequence < 25 bp. Hisat2 was used to align clean reads with the human reference genome Ensemble GRCh38. Expression of mRNAs and lncRNAs was normalized and outputted with StringTie. Then, CIRI2 software was used to predict cir-cRNAs.

Identification of DEmRNAs, DElncRNAs and DEcircRNAs
Ballgown was applied to identify DEmRNAs, DElncRNAs and DEcircRNAs in HS with |log 2 FC| > 1 and P < 0.05. Hierarchical clustering analysis of DEmRNAs, DElncR-NAs and DEcircRNAs was performed with R (https:// www.r-project.org/) package 'pheatmap'. David 6.8 was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for DEmRNAs with P < 0.05.

Discussion
HS, a fibroproliferative disorder, is characterized by excessive deposition of ECM and invasive growth of fibroblasts [12]. In this study, a total of 487 DEmR-NAs, 92 DElncRNAs and 17 DEcircRNAs were identified in HS. In addition, several pathways were identified to be closely associated with HS, including collagen fibril organization, ECM-receptor interaction and PI3K-Akt signaling pathway. According to the    results of functional annotation analysis, ADAMTS14, ACAN and COL11A1 were enriched in collagen fibril organization; COMP was enriched in ECM-receptor interaction; and COMP and COL11A1 were enriched in the PI3K-Akt signaling pathway. Collagen is the main component of interstitial ECM, which is involved in the regulation of various biological processes, such as cell morphology, proliferation, migration, differentiation, apoptosis and carcinogenesis [13]. As a member of minor fibrillar collagens, COL11A1 can be produced by cartilage and a variety of noncartilaginous tissues, including skin [14]. Recent studies have found that COL11A1 is associated with various cancers, such as gastric cancer, ovarian cancer and non-small cell lung cancer [15][16][17]. It has been concluded that COL11A1 expression is a biomarker of human carcinoma-associated stromal cells and carcinoma progression [18]. COL11A1 has been found to be overexpressed in human keloid fibroblasts related to normal skin fibroblasts [19]. Elevated COL11A1 was also observed in scleroderma skin, another condition with extensive fibroblast activation [20]. P4HA3 encodes a component of prolyl 4-hydroxylase, a key enzyme in collagen synthesis. In general, the expression of P4HA3 is very low in normal fetal and adult tissues [21]. In poorly differentiated gastric adenocarcinoma cancer cell line MKN-45 and AGS cells, up-regulated P4HA3 could enhance cell motility and invasiveness [22]. Highly expressed P4HA3 was associated with poor prognosis in gastric cancer [23]. There were no previous reports of the function of P4HA3 in HS. Here, significantly up-regulated COL11A1 and P4HA3 were observed in HS, adding evidence that COL11A1 and P4HA3 may play an important role in HS.
CD1A encodes a member of the CD1 family of transmembrane glycoproteins, which are structurally related to the major histocompatibility complex proteins and form heterodimers with b 2 -microglobulin. The amount of CD1As of the positive dendritic cell was significantly higher in HS than the controlled normal skin [24]. However, CD1A was identified to be the  most significantly down-regulated DEmRNA in this study, indicating that its exact role in HS needs further study to determine. A recent study indicated that LINC01189 was significantly altered in peripheral blood mononuclear cells of patients with rheumatoid arthritis, suggesting that LINC01189 may be a potential biomarker for rheumatoid arthritis [25]. Besides, Yao et al. [26] demonstrated that LINC01189 was down-regulated in hepatitis C virus-infected hepatocellular carcinoma (HCC) tumors and cell lines and may confer a suppression effect on the development of HCC. Other than that, the expression pattern or regulatory effects of LINC01189 in other human diseases have never been elucidated. In other words, LINC01189, the most significantly down-regulated DElncRNA, was first reported to be associated with HS in this study. Aggrecan, encoded by ACAN, is a major proteoglycan component in the ECM of the growth plate and articular cartilage [27]. Mutations in ACAN were reported to be associated with growth defects ranging from mild idiopathic short stature to severe skeletal dysplasias [28]. TNFRSF9, also termed 4-1BB and CD137, is a member of the tumor necrosis factor receptor superfamily, which contributes to the clonal expansion, survival, and development of T cells. It has been suggested that TNFRFS9 expression was a biomarker for tumor-infiltrating lymphocytes in ovarian cancer and melanoma [29]. In addition, TNFRSF9 methylation has been reported to serve as a biomarker in the context of immunotherapies in melanoma [30]. Recently, lncRNA MIR503HG has been suggested to be dysregulated and involved in a variety of human cancers. Qiu et al. [31] suggested that MIR503HG exhibited significant antiproliferation and antimigration/invasion effects on bladder cancer cells. Chuo et al. [32] demonstrated that MIR503HG overexpression inhibits colorectal cancer cell migration and invasion mediated by transforming growth factor-b2. Lin et al. [33] revealed that MIR503HG suppressed nonsmall cell lung cancer progression via negatively regulating Wnt1 expression. ACAN, TNFRSF9 and MIR503HG were identified to be dysregulated in this study, although no previous study linked ACAN, TNFRSF9 and MIR503HG with HS. In addition, ACAN and TNFRSF9 were targets of MIR503HG (MIR503HG/hsa-miR-204-3p/ACAN and MIR503HG/ hsa-miR-431-5p/TNFRSF9) in the ceRNA network, which indicated that MIR503HG may act as a ceRNA to regulate the expression of ACAN and TNFRSF9 in HS.
ADAMTS14, located on chromosome 10q22.1, is a member of the ADAMTS metalloproteinase family, comprised of 19 members, which are known as proteolytic enzymes to catalyze a great variety of substrates in the ECM [34]. The activation of ADAMTS proteinases can exhibit both inhibitory and promotive effects on angiogenesis because the mechanism involved in their regulation of cancer development varies among different members [35]. Sheu et al. [36] implicated the ADAMTS14 gene polymorphism as a predictive factor of HCC. Low cytoplasmic expression of ADAMTS14 has been associated with poor overall survival of patients with oral squamous cell carcinoma, which may be used as a novel biomarker for oral squamous cell carcinoma diagnosis [37]. MEG3, located on chromosome 14q32.3, has been associated with various tumors and regarded as a putative cancer biomarker and treatment target [38]. In the ceRNA network, ADAMTS14 was a target of MEG3 (MEG3/ hsa-miR-6884-5p/ADAMTS14). In the DElncRNA-DEmRNA coexpression network, ADAMTS14 was coexpressed with AC000035.1 (one of the top 10 significantly up-regulated DElncRNAs). Hence we speculated that MEG3 and AC000035.1 may participate in HS via regulating ADAMTS14.
COMP is a fibrillar collagen assembly regulator, which is involved in the assembly and stabilization of the ECM via its interactions with type I and type II collagen and modulates the cellular phenotype during tissue genesis and remodeling [39,40]. Zachou et al. [41] suggested COMP as a biomarker of liver fibrosis in patients with chronic viral hepatitis. Li et al. [39] reported that hepatic stellate cell-derived COMP drives HCC progression by activating mitogenactivated protein kinase kinase 7 (MEK)/mitogenactivated protein kinase (ERK) and PI3K/AKT signaling pathways. Vuga et al. [42] demonstrated that COMP may serve as a biomarker for idiopathic pulmonary fibrosis. Agarwal et al. [43] indicated that COMP is also a constitutive component present in human skin that is deposited by fibroblasts into the ECM of human skin. Agarwal et al. [44] recently demonstrated that COMP deposition is enhanced in the dermis in various fibrotic conditions. ADAM12 encodes a member of the a disintegrin and metalloprotease (ADAM) protein family and is restrictively expressed in normal tissues [45]. It was reported that up-regulated ADAM12 in the central part of keloids may be involved in processes leading to clinical regression [19]. In the DEcircRNA-DEmRNA coexpression network, hsa_circ_0069865, one of the top three downregulated DEcircRNAs that covered the most DEmR-NAs, was coexpressed with COMP and ADAM12, which may suggest hsa_circ_0069865 was involved in HS mediated by COMP and ADAM12.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article.  (D) KEGG pathways. The x axis shows P value of GO terms or KEGG pathways, and the y axis shows GO terms or KEGG pathways. The color scale represented Àlog P value. The x axis shows P value of GO terms or KEGG pathways, and the y axis shows GO terms or KEGG pathways. The color scale represented Àlog P value.