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

Background: Melanoma is the most common and dangerous skin tumor, and its pathogenesis is not fully understood. Methods: To identify the key lncRNAs in melanoma, we reconstructed a global triple network based on the ceRNA theory. GO and KEGG pathway analysis were performed using DAVID. And, we verified our findings through qRT-PCR assay. Results: According to ceRNA theory, 898 DEMs,53 DELs and 16 DEMis were selected to reconstruct the ceRNA network. MALAT1, LINC00943, and LINC00261 were selected as hub gene in the ceRNA network. Conclusions: In this study, we found that MALAT1, LINC00943 and LINC00261 were closely associated with tumorigenesis and development of melanoma, and therefore, could be used as potential diagnostic biomarkers and therapeutic targets.


Background
Melanoma is the most common and dangerous skin cancer. 1,2 Worldwide, cutaneous melanoma occurred in approximately 232,000 (1.7%) cases of all newly diagnosed primary malignant cancers, and resulted in about 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 USA from 1982 to 2011, had shown that the incidence increased about 3% annually, and will further increase at least until 2022. 3 In 2015, there were 3.1 million people with melanoma, resulting in 59,800 deaths. 4 Still, 95,710 cases of melanoma in situ will be newly diagnosed in 2020. 6 The high incidence and high mortality of melanoma make it necessary for us to pay more attention on it. 5,6 Although some achievements have been made in the genetic research of melanoma, more molecules related to diagnosis and treatment have been required to be found.
Tumorigenesis are often resulted from aberrant transcriptomes, including the aberrant levels of coding RNA and non-coding RNA. 7.8,9 It had been proved that lncRNA have a variety of effects, including regulation of gene transcription, regulation of post-transcription and epigenetic regulation and more. [10][11][12][13] In addition, dysregulation of lncRNA had been observed in various diseases. [14][15][16][17][18] Unfortunately, the function of lncRNA are more difficult to identified in comparison with those of coding RNAs. Until now, only few lncRNA was identified as crucial factor in tumorigenesis and development of melanoma, including ZNNT1, THOR and SAMMSON. [19][20][21] Thus, how to locate them and define its functions are the challenge of current research.
The effect of miRNA on malignancies had been well verified in many ways. Studies had suggested that lncRNAs can regulate miRNAs abundance by binding and sequestering them. 22 Thus, it gave us a hint to study the function of lncRNAs, which is, by studying the interaction among lncRNAs, mRNAs and miRNAs. In 2011, the competitive endogenous RNA (ceRNA) hypothesis, a novel regulatory mechanism between non-coding RNA and coding RNA had been proposed. [23][24][25] This theory indicated that any RNA transcript harboring miRNA-response elements (MREs) can sequester miRNAs from other targets sharing the same MREs, and thereby regulate their expressions. [23][24][25] That is to say, the RNA transcripts which can be cross regulated by each other can be biologically predicted according to their common MREs. 24.26 Evidences had been shown that ceRNAs existed in several species and contexts, and might play an important role in various biological processes, such as tumorigenesis. 25 Systematic analysis of ceRNA network had been well performed in multiple tumors, such as gastric cancer, bladder cancer, and ovarian cancer, contributing to a better understanding of tumorigenesis and facilitate the development of lncRNA-directed diagnostics and therapeutics against this disease. [27][28][29] Unfortunately, however, such functional interactions have not yet been elucidated in melanoma.
In this study, ceRNA networks had been reconstructed, in order to reveal the potential diagnosis and therapeutic molecules of melanoma.

Raw data
Human melanoma miRNA expression data were download from NCBI GEO database (GEO (http://www.ncbi.nlm.nih.gov/geo), 32 including GSE24996, GSE35579, GSE62372, which were array-based datasets. The GSE24996 dataset consists of 8 benign nevus tissue samples and 15 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 download from NCBI GEO database (GSE112509), which is sequence-based dataset. The GSE112509 dataset consist of 23 benign nevus tissue samples and 57 primary melanoma tissue samples, respectively.

Prediction of target lncRNAs and mRNAs
To predict the target lncRNAs and mRNAs through DEMis, starbase (starbase.sysu.edu.cn) was used in our study. 36 Multiple lncRNA/mRNA-predicting programs (PITA, RNA22, miRmap, DIANA-microT, miRanda, PicTar and TargetScan) were used in starbase. 36 In order to make our results more accurate, only when the target mRNA was predicted in at least four predicted programs on starbase, it would be chosen as the predicted target mRNA. Then, these predicted target lncRNAs and mRNAs were merged with DEMs and DELs, respectively.

Reconstruction of ceRNA network
The ceRNA network have been reconstructed based on ceRNA theory, and as follows 24 : (1) Expression correlation between DELs and DEMs was evaluated using the Pearson correlation coefficient (PCC). The DELs-DEMs pairs with PCC > 0.4 and P-value < 0.01 were considered as co-expressed lncRNA-mRNA pairs. (2) Both lncRNAs and mRNAs in the pairs were negatively correlated with their common miRNAs. To reconstruct our key ceRNA sub-network, we first select hub gene according to the node degrees of 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 neighbors, by using DAVID (version 6.8, https://david.ncifcrf.gov/) 30,31 .

Selection of melanoma sample
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 World Medical Association (Declaration of Helsinki). Eligible patients were from The First Affiliated Hospital, Sun yat-sen University (Guangzhou, Guangdong, China) or Cancer Center of Guangzhou Medical University (Guangzhou, Guangdong, China). The eligible patients for study had to meet the following criteria: (1) histologically confirmed as melanoma; (2) received no radiotherapy, chemotherapy or biotherapy before surgery. Exclusion criteria include: (1) with previous malignancies; (2) with concomitant malignancies; (3) with serious active infection; (4) experiencing pregnancy or lactation. Melanoma tissues and its adjacent normal tissues, were collected from each patient. The adjacent normal tissue was taken more than 2cm away from the edge of the tumor, which was confirmed as normal skin tissue by at least three pathologists.

Selection of healthy sample
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 World Medical Association (Declaration of Helsinki). Eligible patients were from The First Affiliated Hospital, Sun yat-sen University (Guangzhou, Guangdong, China). Healthy patients between 18 and 50 years of age who were scheduled to undergo split-thickness skin grafting were included in the study. Exclusion criteria include: (1) with previous malignancies; (2) with concomitant malignancies; (3) with serious active infection; (4) experiencing pregnancy or lactation; (5) with diabetes or other skin diseases.

RNA isolation and qRT-PCR
The total RNA was extracted from all samples by the use of Trizol reagent (Invitrogen, USA). The OD value (260/280) of all RNA extracted samples was greater than 1.8. For each replicate, complementary DNA (cDNA) was synthesized from 2μg RNA using GoScript Reverse Transcription System (Promega, USA). The qRT-PCR reaction comprised 10μl of GoTaq qPCR Master Mix (2×) (Promega, USA), 2μl of diluted cDNA template (1:10) and 10μM of each primer contributing to a total volume of 20μl. Reactions were run in ABI 7500 real-time PCR system (Applied Biosystems, USA) under the following conditions: 95°C for 10 mins, and 40 cycles of 95°C for 15s and 60°C for 60s. Melting curves were derived for every reaction to insure a single product. Relative gene expression was evaluated according to ddCT method, using human GAPDH gene as endogenous controls for RNA load and gene expression in analysis. All experiments were performed in triplicate. The GraphPad Prism 8 (GraphPad Software, USA) was used to output figures.

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

Hub gene selection
According to the node degree in ceRNA network, we found that three lncRNAs, MALAT1, LINC00943, 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 3. Therefore, these three lncRNAs might play an essential role in melanomagenesis, which might be considered as the key lncRNAs. The MALAT1LINC00943 and LINC00261 sub-network show that GO-BP analysis revealed 75, 69 and 42 enriched terms respectively, which are also significantly enriched in regulation of transcription such as positive regulation of transcription from RNA polymerase II promoter, positive regulation of transcription (DNA-templated), transcription from RNA polymerase II promoter and more. Pathways analysis reveal that 19, 13 and 7 pathways are significantly enriched in tumor-related pathways, including pathway in cancer, focal adhesion, PI3K-Akt signaling pathway and more.

Expression of MALAT1LINC00943 and LINC00261 are higher in tumor tissues
To confirm the expression of MALAT1LINC00943 and LINC00261 in melanoma tissues, we had evaluated MALAT1LINC00943 and LINC00261 expression level in cancer tissues from 12 melanoma patients and 3 healthy samples via qRT-PCR, as shown in Fig. 7. The results show that the expression of MALAT1LINC00943 and LINC00261 are significantly higher in tumor tissues compared with healthy samples (p = 0.0243, p = 0.0005, p < 0.0001, respectively). Also, the expression of MALAT1LINC00943 and LINC00261 are significantly higher in tumor tissues compared with its adjacent normal tissue (p = 0.0002, p < 0.0001, p < 0.0001, respectively). However, no significant difference was observed between healthy samples and adjacent normal skin tissues in expression of MALAT1LINC00943 and LINC00261 (p = 0.366, p = 0.379, p = 0.262, respectively). The results are consistent with that discussed above. Thus, the expression of MALAT1LINC00943 and LINC00261 are increased in melanoma and may responsible for the tumorigenesis of melanoma.

Discussion
Melanoma is the most common skin cancer. 1,2 The high incidence and high mortality of melanoma make it necessary for us to pay more attention on it. 5,6 Over the past few years, great efforts have been made to explore the molecular mechanisms of it. 6 Previous studies have been focused on the protein-coding genes and miRNAs, however, it rarely involves in lncRNAs. In this study, three lncRNAs, MALAT1, LINC00943 and LINC00261, were identified according to the reconstructed ceRNA network.
With the increasing attention in the role of lncRNA, lncRNA have shown superior potential over protein-coding genes as a biomarker for diagnosis, prognosis and treatment. [39][40][41][42] Among these key lncRNAs found in this study, MALAT1 has been demonstrated to be related to various malignant tumors. [43][44][45][46][47] Studies had confirmed that MALAT1 is a valuable prognostic marker and a promising therapeutic target in lung cancer metastasis. 43, 44 A study had also suggested that MALAT1 play an important role in tumor progression and could be served as a promising therapeutic target. 45 Through the study of whole-genome mutational landscape and characterization of noncoding and structural mutations in liver cancer, Fujimoto A. and colleges discovered MALAT1 closely related to liver carcinogenesis. 46 In addition, a study revealed a novel mechanism of MALAT1-regulated autophagy-related chemoresistance in gastric cancer. 47 At present, it is believed that MALAT1 is mainly responsible for regulating the proliferation, migration and invasion of tumor cells. According to our findings, which indicated that MALAT1 might also be a crucial factor in the tumorigenesis and development of melanoma. In this sub-network, we found that there are nine lncRNA-miRNA pairs, including miRNA-378a-3p, miRNA-23b-3p, miRNA-224-5p, miRNA-204-5p, miRNA-205-5p, miRNA-200c-3p, miRNA-200b-3p, miRNA-149-5p, miRNA-211-5p. Among them, study have demonstrated that MALAT1 can regulate chemoresistance via miRNA-23b-3p sequestration in gastric cancer. 47 In ovarian cancer, a study suggested that MALAT1-miRNA-211-5p may act as a key mediator in prevention of this disease. 48 MALAT1 also involved in promoting renal cell carcinoma through interaction with miRNA-205-5p. 49 Studies have confirmed that MALAT1 functions in liver and lung cancer through miRNA-204-5p. 50,51 In addition, targeting MALAT1/miRNA-200c-3p axis in xenograft endometrial carcinoma model can greatly inhibit tumor growth. 52 Moreover, studies had been illustrated that these miRNAs were closely related to melanoma in serval ways.
miRNA-378a-3p can regulate oncogene PARVA expression in melanoma, preventing its progression. 53 miRNA-23b-3p was shown as a tumor suppressor gene in melanoma. 54 miRNA-224-5p can be regulated by E2F1 to drive EMT through TXNIP downregulation in melanoma, and it can inhibit uveal melanoma cell proliferation, migration, invasion by targeting PIK3R3/AKT3. 55,56 miRNA-204-5p, known as tumor suppressor gene in melanoma, was associated with CDKN2A pathway and NRAS gene, and was contributed to BRAF inhibitor resistance. 54.57,58 miRNA-205-5p suppresses proliferation and induces senescence via regulation of E2F1 in melanoma. 54.59-61 miRNA-200b/c-3p act as potential diagnostic and prognostic markers for melanoma. 62-64. Upregulation of miRNA-149-5p, directedly regulated by p53, results in increased expression of Mcl-1 and resistance to apoptosis in melanoma cells. 65 Most importantly, studies have confirmed that miRNA-211-5p play a major role as tumor suppressor via various targets in melanoma. 54,58,62,66,67 As these miRNAs we predicted were closely associated with melanoma, we believed that MALAT1 may play a crucial role in tumorigenesis and development of melanoma through competitive interactions with these miRNAs, and subsequently alter expression of downstream mRNAs.
LINC00261 is also significant based on ceRNA network, but its function remains controversial until now. We found that 43 GO-BP terms and 7 KEGG pathways were enriched, based on the sub-network. One of these pathways, PI3K/Akt signaling pathway, is proved to play a critical role in tumorigenesis, 69 especially in melanoma. 70 Also, study have demonstrated that LINC00261 promotes cancer cell proliferation and metastasis in human choriocarcinoma. 71 However, LINC00261 have shown a great capacity in improving the chemotherapeutic response and survival of patients with esophageal cancer. 72 And, in gastric cancer, LINC00261 can suppress the tumor metastasis by regulating epithelial-mesenchymal transition (EMT). 73 Moreover, LINC00261 can block cellular proliferation by activating the DNA damage response. 74 LINC00261 may affect the biological behavior of different tumors in different ways. Therefore, it is very essential to further explore the role of LINC00261 in different tumors. On the other hand, five miRNAs, including miRNA-23b-3p, miRNA-211-5p, miRNA-205-5p, miRNA-140-3p and miRNA-125b-5p were interacted with LINC00261 according to LINC00261-miRNA-mRNA subnetwork. Similarly, no connection between LINC00261 and these miRNAs have been discovered yet, but these miRNAs were also demonstrated to be associated with melanoma. The role of miRNA-23b-3p, miRNA-211-5p, miRNA-205-5p, and miRNA-125b-5p in melanoma were discussed above. miRNA-140-3p was reported to be regulated by MALAT1 in uveal melanoma cell, 75 while its role in cutaneous melanoma is still unknown.
Three of the 16 predicted miRNAs were not associated with MALAT1, LINC00943 and LINC00261, including miRNA-21-5p, miRNA-20b-5p and miRNA-424-5p. They closely related to SGMS1.AS1, EPB41L4A.AS1 and SNHG1, on the other hand. Little was known about miRNA-424-5p in melanoma, while studies have suggested that miRNA-20b-5p may inhibit tumor metastasis via regulation of PAR-1 receptor in melanoma cells, 76 and miRNA-21 may regulate melanoma cell proliferation, migration, apoptosis through ERK/NF-κB signaling pathway by targeting SPRY1, PDCD4 and PTEN. 77,78 All in all, we reconstructed a ceRNA network, which for the first time enables an overall view and analysis of the lncRNA-associated ceRNA mediated genes in the tumorigenesis and development of melanoma at a system-wide level based on ceRNA theory. Our study discovered that serval lncRNAs, including MALAT1, LINC00943 and LINC00261, might play an essential role in skin cutaneous melanoma. Finally, we verified our findings through q(RT)-PCR assay. This study will advance our understanding of the tumorigenesis and development of melanoma from the perspective of lncRNAs and provide some novel lncRNAs as candidate diagnostic biomarkers or potential therapeutic targets for melanoma. Further studies are required to explore the biological functions and molecular mechanisms of MALAT1, LINC00943 and LINC00261 in melanoma.

Conclusions
This study advances our understanding of tumorigenesis and development in melanoma from the perspective of lncRNAs, and provides three novel lncRNAs as candidate diagnostic biomarkers or potential therapeutic targets for melanoma. Further studies are required to verify the role of MALAT1, LINC00943 and LINC00261 in melanoma.

Declarations
Ethics approval and consent to participate: 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 World Medical Association (Declaration of Helsinki).

Tables
Due to technical limitations, all tables are available as downloads in the Supplementary Files section.     The expression level of MALAT1 (A), LINC00943 (B) and LINC00261 (C) in normal skin, adjacent normal skin and melanoma tissues.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.