Reconstruction of lncRNA-miRNA-mRNA network based on competitive endogenous RNA reveals functional lncRNAs in skin cutaneous melanoma

Human skin cutaneous melanoma is the most common and dangerous skin tumour, but its pathogenesis is still unclear. Although some progress has been made in genetic research, no molecular indicators related to the treatment and prognosis of melanoma have been found. In various diseases, dysregulation of lncRNA is common, but its role has not been fully elucidated. In recent years, the birth of the “competitive endogenous RNA” theory has promoted our understanding of lncRNAs. To identify the key lncRNAs in melanoma, we reconstructed a global triple network based on the “competitive endogenous RNA” theory. Gene Ontology and KEGG pathway analysis were performed using DAVID (Database for Annotation, Visualization, and Integration Discovery). Our findings were validated through qRT-PCR assays. Moreover, to determine whether the identified hub gene signature is capable of predicting the survival of cutaneous melanoma patients, a multivariate Cox regression model was performed. According to the “competitive endogenous RNA” theory, 898 differentially expressed mRNAs, 53 differentially expressed lncRNAs and 16 differentially expressed miRNAs were selected to reconstruct the competitive endogenous RNA network. MALAT1, LINC00943, and LINC00261 were selected as hub genes and are responsible for the tumorigenesis and prognosis of cutaneous melanoma. MALAT1, LINC00943, and LINC00261 may be closely related to tumorigenesis in cutaneous melanoma. In addition, MALAT1 and LINC00943 may be independent risk factors for the prognosis of patients with this condition and might become predictive molecules for the long-term treatment of melanoma and potential therapeutic targets.


Background
Human skin cutaneous melanoma (SKCM) is the most common and dangerous type of skin tumour [1,2]. Worldwide, approximately 232,000 (1.7%) cases of cutaneous melanoma are reported among all newly diagnosed primary malignant cancers, and this disease results in approximately 55,500 cancer deaths (0.7% of all cancer deaths) [1,3]. The incidence of melanoma in Australia, New Zealand, Norway, Sweden, the UK, and the USA from 1982 to 2011 has shown increases of approximately 3% annually and will further increase until 2022 [3]. In 2015, there were 3.1 million people with melanoma, resulting in 59,800 deaths [4]. Nevertheless, 95,710 cases of melanoma in situ will be newly diagnosed in 2020 [5]. The high incidence and high mortality of melanoma indicate that researchers must further study this disease. Although some achievements have been made in the genetic research of melanoma, markers related to diagnosis and treatment are needed.
Tumorigenesis often results from aberrant transcriptomes, including aberrant levels of coding RNA and noncoding RNA [6][7][8]. It has been proven that lncRNAs have various effects, including regulation of gene transcription, post-transcriptional regulation and epigenetic regulation [9][10][11][12]. In addition, dysregulation of lncRNAs has been observed in various diseases [13][14][15][16]. Unfortunately, the functions of lncRNAs are more difficult to identify than those of coding RNAs. Until now, only a few lncRNAs have been identified as crucial factors in the tumorigenesis and development of melanoma, including ZNNT1, THOR and SAMMSON [14,15,17]. Thus, how to locate them and define their functions is a challenge of current research.
The effect of miRNAs on malignancies has been verified in many ways. Studies have suggested that lncRNAs can regulate miRNA abundance by binding and sequestering them [18]. Thus, we aimed to study the function of lncRNAs by studying the interactions among lncRNAs, mRNAs and miRNAs. In 2011, the competitive endogenous RNA (ceRNA) hypothesis proposed a novel regulatory mechanism between noncoding RNA and coding RNA [19][20][21]. This theory indicated that any RNA transcript harbouring miRNA-response elements (MREs) can sequester miRNAs from other targets sharing the same MREs and thereby regulate their expression [19][20][21]. That is, the RNA transcripts that can be cross regulated by each other can be biologically predicted according to their common MREs [20,22]. Evidence has shown that ceRNAs exist in several species and contexts and might play an important role in various biological processes, such as tumorigenesis [21]. Systematic analysis of the ceRNA network has been performed in multiple tumours, such as gastric cancer, bladder cancer, and ovarian cancer, contributing to a better understanding of tumorigenesis and facilitating the development of lncRNA-directed diagnostics and therapeutics against this disease [23][24][25]. Unfortunately, however, such functional interactions have not yet been elucidated in melanoma.
In this study, we used bioinformatics methods to construct the ceRNA network of cutaneous melanoma and to identify the key lncRNAs involved in melanomagenesis. Through the reconstruction of a ceRNA network, we identified and verified that the key ceRNA molecules play a crucial role in the tumorigenesis and prognosis of SKCM. (Work flow was shown in Fig. 1).

Raw data
Human melanoma miRNA expression data were downloaded from the NCBI GEO database (GEO (http://   www.ncbi.nlm.nih.gov/geo) [26], including GSE24996, GSE35579, and GSE62372, which are array-based datasets. The GSE24996 dataset consists of 8 benign nevus tissue samples and 23 primary melanoma tissue samples. The GSE35579 dataset consists of 11 benign nevus tissue samples and 20 primary melanoma tissue samples. The GSE62372 dataset consists of 9 benign nevus tissue samples and 92 primary melanoma tissue samples. mRNA and lncRNA expression data were also downloaded from the NCBI GEO database (GSE112509), which is a sequence-based dataset. Venn diagram: (a) DEMis were selected with |log2FC| > 1 and adjusted P-value < 0.05 among the non-coding RNA profiling sets, GSE24996, GSE35579 and GSE62372. The candidates 18 miRNAs were shared in at least two datasets. b DEMs were selected by intersecting mRNAs predicted by DEMis through starbase and differential expressed mRNAs in GSE112509. c DELs were selected by intersecting lncRNAs predicted by DEMis through starbase and differential expressed lncRNAs in GSE112509. (These images were produced by R version 3.4.2)  Table 3. (Fig. 4a was produced by Cytoscape version 3.7.1)  The GSE112509 dataset consists of 23 benign nevus tissue samples and 57 primary melanoma tissue samples.

Prediction of target lncRNAs and mRNAs
For prediction of the target lncRNAs and mRNAs through DEMis, starBase (starbase.sysu.edu.cn) was used in our study [30]. Multiple lncRNA/mRNApredicting programmes (PITA, RNA22, miRmap, DIANA-microT, miRanda, PicTar and TargetScan) were used in starBase [30]. For accuracy, only when the target mRNA was predicted in at least four predicted programmes on starBase would it be chosen as the predicted target mRNA. Then, these predicted target lncRNAs and mRNAs were merged with DEMs and DELs, respectively.

Reconstruction of the ceRNA network
The ceRNA network was reconstructed based on ceRNA theory [20] and as follows: (1) Expression correlation between DELs and DEMs was evaluated using the Pearson correlation coefficient (PCC). The DEL-DEM pairs with PCC > 0.4 and P-value < 0.01 were considered coexpressed lncRNA-mRNA pairs. (2) Both lncRNAs and mRNAs in the pairs were negatively correlated with their common miRNAs. (3) The ceRNA network was reconstructed and visualized using Cytoscape (version 3.7.1, https://cytoscape.org/) [31,32].

Hub gene selection and reconstruction of key ceRNA subnetworks
To reconstruct our key ceRNA subnetwork, we first selected hub genes according to the node degrees of the ceRNA network we reconstructed above by calculating the number of lncRNA-miRNA and miRNA-mRNA pairs. For these key lncRNAs, GO-BP, GO-CC, GO-MF and KEGG pathway annotation were performed according to their first mRNA neighbours by using DAVID (version 6.8, https://david.ncifcrf.gov/) [33,34].

Sample selection for qRT-PCR validation
To validate findings in the ceRNA network, we selected the top three hub genes to determine their expression in cutaneous melanoma and skin tissues. Twelve patients with cutaneous melanoma and three healthy patients were included in this study. The study protocol was approved by the Ethics Committee of The First Affiliated Hospital, Sun Yat-sen University. All patients provided written informed consent in compliance with the code of ethics of the World Medical Association (Declaration of Helsinki). The eligible patients for this study had to meet the following criteria: (1) histologically confirmed as melanoma; (2) received no radiotherapy, chemotherapy or biotherapy before surgery. The exclusion criteria were as follows: (1) previous malignancies; (2) concomitant malignancies; (3) serious active infection; and (4) pregnancy or lactation. Eligible cutaneous melanoma patients were from The First Affiliated Hospital, Sun Yat-sen University (Guangzhou, Guangdong, China) or the Cancer Center of Guangzhou Medical University (Guangzhou, Guangdong, China). Each tumour sample was matched with adjacent apparently normal tissues removed during the same operation. Frozen sections were made from these tissues and examined by at least three pathologists. The clinicopathological features of twelve skin cutaneous melanoma patients (51.67 ± 14.57 years old) for qRT-PCR validation are shown in Table 1. Three healthy patients from The First Affiliated Hospital, Sun Yat-sen University (Guangzhou, Guangdong, China) were included in this study. These patients were scheduled to undergo split-thickness skin grafting due to deep partial burn wounds. Each normal skin sample was obtained from the donor site. All the samples were frozen immediately after the operation and were stored in liquid nitrogen until RNA isolation.   Table 5. (Fig. 5a

KEGG pathway and GO enrichment analysis of lncRNAs based on the ceRNA network
We used DAVID to analyse the biological classification of DEMs according to method 2.5. The results of the top 15 significant GO terms and KEGG pathways are shown in Table 3 and Fig. 4b-e. Sixty pathways were significantly enriched through KEGG pathway analysis, including the PI3K-Akt signalling pathway, focal adhesion, proteoglycans in cancer, pathway in cancer and, most importantly, melanomagenesis. The results of GO-BP analysis revealed 172 enriched terms, particularly in the regulation of transcription, such as positive regulation of transcription from the RNA polymerase II promoter, positive regulation of transcription (DNA-templated), and transcription from the RNA polymerase II promoter.

Hub gene selection
According to the node degree in the ceRNA network, we found that three lncRNAs, MALAT1, LINC00943, and LINC00261, had the highest number of lncRNA-miRNA and miRNA-mRNA pairs, suggesting that these three lncRNAs could be chosen as hub nodes, and the results are shown in Table 4. Therefore, these three lncRNAs might play an essential role in melanomagenesis and might be considered key lncRNAs.     Fig. 5b-e, Fig. 6b-e, Fig. 7b-e, and Tables 5, 6, 7.

Expression of MALAT1, LINC00943 and LINC00261 is higher in tumour tissues
To confirm the expression of MALAT1, LINC00943 and LINC00261 in melanoma tissues, we evaluated the MALAT1, LINC00943 and LINC00261 expression levels in the cancer tissues from 12 melanoma patients (see Table 1) and 3 healthy tissues via qRT-PCR, as shown in Fig. 8. The results showed that the expression of MALAT1, LINC00943 and LINC00261 was significantly higher in the tumour tissues than in the healthy tissues (p = 0.0243, p = 0.0005, p < 0.0001, respectively). Additionally, the expression of MALAT1, LINC00943 and LINC00261 was significantly higher in the tumour tissues than in the adjacent normal tissues (p = 0.0002, p < 0.0001, p < 0.0001, respectively). However, no significant difference was observed between the healthy tissues and the adjacent normal skin tissues in the expression of MALAT1, LINC00943 and LINC00261 (p = 0.366, p = 0.379, p = 0.262, respectively). The results are consistent with those discussed above. Thus, the expression of MALAT1, LINC00943 and LINC00261 is increased in melanoma and may be responsible for the tumorigenesis of melanoma.
MALAT1 and LINC00943 are independent risk factors for the prognosis of cutaneous melanoma A univariate Cox regression model for survival analysis of age, sex and stage was performed, and the results are shown in Supplementary Table 3. Then, the multivariate Cox regression model for survival analysis of MALAT1, LINC00943, and LINC00261 was performed. The results showed that the overall survival time and disease-free survival time of the patients with MALAT1 or LINC00943 CNV deficiency were significantly lower than those without it, and the difference was significant (details are shown in Table 8 and Fig. 9a-d), suggesting that MALAT1 and LINC00943 are independent risk factors for the prognosis of cutaneous melanoma. Although the overall survival time and disease-free survival time of patients with LINC00261 deletion were lower than those without it, the difference was not significant (p = 0.535, p = 0.694) (details are shown in Table 8 and Fig. 9e-f).

Discussion
In this study, three lncRNAs, MALAT1, LINC00943 and LINC00261, were identified according to the reconstructed ceRNA network. Among these key lncRNAs found in this study, MALAT1 has been demonstrated to be related to various malignant tumours [40][41][42][43][44]. Studies have confirmed that MALA T1 is a valuable prognostic marker and a promising therapeutic target in lung cancer metastasis [40,41]. A study also suggested that MALAT1 plays an important role in tumour progression and could serve as a promising therapeutic target [42]. Through the study of the whole-genome mutational landscape and characterization of noncoding and structural mutations in liver cancer, Fujimoto A. and colleagues discovered that MALAT1 is closely related to liver carcinogenesis. 46 In addition, a study revealed a novel mechanism of MALAT1-regulated autophagy-related chemoresistance in gastric cancer [44]. At present, it is believed that MALAT1 is mainly responsible for regulating the proliferation, migration and invasion of tumour cells. According to our findings, MALAT1 might also be a crucial factor in the tumorigenesis and development of melanoma.
Seven KEGG pathways were enriched based on the LINC00261 subnetwork. One of these pathways, the PI3K/Akt signalling pathway, has been proven to play a critical role in tumorigenesis [66], especially in melanoma [67]. Additionally, a study has demonstrated that LINC00261 promotes cancer cell proliferation and metastasis in human choriocarcinoma [68]. However, LINC00261 has shown a strong capacity in improving the chemotherapeutic response and survival of patients with oesophageal cancer [69]. In gastric cancer, LINC00261 can suppress tumour metastasis by regulating epithelialmesenchymal transition [70]. Moreover, LINC00261 can block cellular proliferation by activating the DNA damage response [71]. LINC00261 may affect the biological behaviour of different tumours in different ways. Therefore, it is essential to further explore the role of LINC00261 in different tumours. However, five miRNAs, miRNA-23b-3p, miRNA-211-5p, miRNA-205-5p, miRNA-140-3p and miRNA-125b-5p, interacted with LINC00261 according to the LINC00261-miRNA-mRNA subnetwork. Similarly, no connection between LINC00261 and these miR-NAs has been discovered yet. The roles of miRNA-23b-3p, miRNA-211-5p, miRNA-205-5p, and miRNA-125b-5p in melanoma are discussed above. miRNA-140-3p was reported to be regulated by MALAT1 in uveal melanoma cells [72]. The multivariate Cox regression model for survival suggested that LINC00261 was not a risk factor for the prognosis of SKCM, however, the median overall survival and disease-free survival time for patients with LINC00261 CNV deficiency were significantly lower than those without LINC00261 CNV deficiency (17.03 m vs 61.05 m, 13.50 vs 25.02).

Conclusions
This study advances our understanding of tumorigenesis and development in cutaneous melanoma from the perspective of the ceRNA theory. In addition, MALAT1 and LINC00943 may be independent risk factors for the prognosis of patients with cutaneous melanoma and might become predictive molecules for the long-term treatment of melanoma and potential therapeutic targets. Further studies are required to validate the role of MALAT1, LINC00943 and LINC00261 in cutaneous melanoma.