Identification of Long Non-coding and Messenger RNAs Differentially Expressed Between Primary and Metastatic Melanoma

Purpose: Melanoma is the most aggressive and life-threatening cutaneous cancer. To explore new treatment strategies, it is essential to identify the mechanisms underlying melanoma tumorigenesis and metastasis.

Methods: In the current study, we demonstrated altered expression of long non-coding RNA (lncRNA) and messenger RNA (mRNA) in melanoma using data from the Cancer Genome Atlas (TCGA) database. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and protein–protein interaction (PPI) analyses were conducted. We also constructed a functional lncRNA-mRNA regulatory network and Kaplan-Meier analysis.

Results: We identified 246 differentially expressed (DE) lncRNAs and 856 DEmRNAs. A total of 184 DElncRNAs and 428 DEmRNAs were upregulated in metastatic melanoma, while all others were downregulated. Additionally, we investigated the co-expression pattern of 363 genes, among which 26 upregulated lncRNAs, 9 down- regulated lncRNAs, 49 upregulated mRNAs and 151 downregulated mRNAs were identified as being co-expressed with others. Survival analysis suggested high levels of 14 lncRNAs and 10 mRNAs may significantly increase or decrease overall survival. These differentially expressed genes are also potentially prognostic in melanoma.

Conclusion: Our findings observe potential roles for lncRNAs and mRNAs during melanoma progression and provide candidate biomarkers for further studies.

Introduction

Cutaneous melanoma is the most aggressive type of skin cancer, and melanoma metastasis accounts for the vast majority of patient deaths from skin cancer. For early stage melanoma, surgical excision is used as the standard approach to remove localized lesions and results in a good prognosis. However, the prognosis of patients with metastatic melanoma is extremely poor. The 5-year survival rate of metastatic melanoma is only 5–10%, while that of non-metastatic melanoma is estimated to be approximately 90%. Although immunotherapy and molecularly targeted drugs are emerging as attractive new therapeutic approaches (Reddy et al., 2016), their utility is restricted to only a portion of patients who often acquire drug resistance, and the therapeutic options against late-stage melanoma remain very limited. A comprehensive understanding of the genetic and epigenetic underpinnings of melanoma initiation and progression would facilitate the identification of new biomarkers and therapeutic targets for metastatic melanoma.

Recently, actively transcribed lncRNAs have been demonstrated to be involved in complex cancer genome regulatory networks rather than representing merely extraneous transcription as originally believed. lncRNAs are defined as endogenous cellular RNA transcripts longer than 200 nucleotides without protein-coding capacity. Their structures resemble mRNAs, but they exhibit lower levels of expression and more tissue-specific expression patterns than protein-coding genes (Cabili et al., 2011). lncRNAs play extensive functions encompassing biological functions including biological evolution, cell proliferation, morphogenesis, cell cycle, immune and inflammatory response, and restrainting miRNA, mRNA, and proteins from combining with their intended targets (Cathcart et al., 2015; Chen et al., 2017; Yu et al., 2018). The expression of lncRNAs varies significantly (Yan et al., 2015) in different tissues and stages of development.

To date, the roles of most lncRNAs have not been elucidated in spite of their relatively huge representation in genomes. It is well known that co-expressed genes are more likely to be co-regulated and functionally related (Stuart et al., 2003). Thus, identifying co-expressed protein-coding genes may reveal functions to associated lncRNAs through use of publicly available genome-wide datasets and computational methods (Signal et al., 2016).

In recent years, TCGA¹ has generated comprehensive, multi-dimensional maps of the crucial genomic changes in numerous cancers, which help to provide an understanding of how such changes interact to drive diseases. Weighted correlation network analysis (WGCNA), a powerful correlation network method, is used to reconstruct gene co-expression network relationships between differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEmRNAs). Gene Ontology Enrichment Analysis Software Toolkit (GOEAST), PANTHER, and the database for annotation, visualization, and integrated discovery (DAVID) (Huang da et al., 2009) combined with the Search Tool for the Retrieval of Interacting Genes (STRING) database² are commonly used to explore the biological functions of the genomic changes. This method has been successfully applied in various biological contexts to identify candidate biomarkers or therapeutic targets, e.g., in cancer, mouse genetics, yeast genetics, and the analysis of brain imaging data. In the present study, we used a computational strategy to perform a systematic study of lncRNAs significantly altered in the transition from primary to metastatic melanoma, and constructed a functional lncRNA-mRNA regulatory network aiming to contribute to a comprehensive understanding of the differentially expressed genes (DEGs) in melanoma and their function. We identified novel melanoma-associated lncRNAs and predicted their potential pathological and biological roles in cancer development and gene regulation. This work may provide some powerful evidence to guide subsequent experimental studies on the altered genes in melanoma as candidate targets for treatment development or as biomarkers.

Materials and Methods

Differentially Expressed lncRNA and mRNA Analyses

Transcriptome sequencing data by RNA-seq were downloaded from the publically available TCGA database, including 103 primary melanoma frozen tissue and 368 metastatic melanoma tissue datasets. All samples were collected from individuals of different races and genders. This dataset comprised of called gene counts for 57,073 mRNAs and non-coding (nc)RNAs. To obtain genome-wide lncRNA and mRNA expression profiles, normalized expression data were subsequently analyzed for DElncRNAs and mRNAs using Bioconductor package (limma, version 3.4.0) in R software (version 3.4.0) with default parameters (Ritchie et al., 2015). A total of 246 DElncRNAs and 856 DEmRNAs remained in the differential expression analysis using “edgeR” in R. Statistical significance was defined as |log₂foldchange| >2 and P-value < 0.01. The FDR method by the Benjamini and Hochberg method was applied to correct for multiple testings. Hierarchical clustering, shown in the format of a heatmap and a volcano plot, of the DElncRNA and DEmRNA expression profiles was performed by using ”RColorBrewer,” “pheatmap,” and “plot” in R language.

Co-expression Network Construction of the lncRNA-mRNA Regulatory Network

A DElncRNA-DEmRNA co-expression network was generated to explore the key roles of DElncRNAs and DEmRNAs in the progression of melanoma. As a systems biology method, gene co-expression network analysis was performed by the WGCNA package (Langfelder and Horvath, 2008) to construct a co-expression network between the onco-lncRNAs and their target mRNAs, based on the signed Pearson Correlation Coefficient between their normalized expression levels as provided by Cuffnorm (Trapnell et al., 2012). As to achieve a scale-free topology, we set β = 7 and modified the pairwise correlation into an adjacency matrix of connection strengths through a soft-thresholding approach. Only edges with weight above a threshold of 0.8 are displayed. In the co-expression network, nodes represented DEGs, and the correlation of gene expression pattern was defined as connectivity degree among genes (Ghazalpour et al., 2006). The network was visualized using Cytoscape 3.5.1.

Bioinformatics and PPI Network Analysis

Gene Ontology analysis is a functional analysis for exploring DEmRNAs with GO categories. The predicted target genes above were input into David³, GOEAST⁴ and PANTHER⁵, which utilized GO to identify coexisting biological process, cellular component, and molecular function represented in the gene profile. Furthermore, we also used the DAVID database to analyze the potential functions of these genes in the pathways (P < 0.05). The results were in accordance with those from KEGG⁶.

Protein–protein interaction information of DEmRNAs was acquired from the STRING database⁷. Edges with weight above a threshold of 0.4 are displayed. Hub genes were found by using “barplot” in R language.

Survival Rate Analysis

Furthermore, we investigated whether the co-expressed DElncRNAs-DEmRNAs were correlated with potential good or poor overall survival (OS). Clinical data from TCGA with survival time were collected. The gene expression was labeled as high or low using a dichotomy method (Tian et al., 2017). OS rate analysis was performed by using “hash” and “survival” in R to examine the association between the gene expression levels and the overall lifetime (Gentleman et al., 2006). Differences in the OS between the two groups were estimated and compared by the Kaplan–Meier method with a log-rank test. P-values less than 0.05 were considered significantly different.

Receiver Operating Characteristic (ROC) Analyses

In order to evaluate the diagnostic value of candidate DElncRNAs, ROC analyses were performed using the “survivalROC” package in R language. AUC under a binomial exact confidence interval was calculated to generate the ROC curve. Univariate Cox and multivariate Cox analyses were evaluated separately.

Results

Identification of DElncRNAs and DEmRNAs

A gene expression dataset generated by RNA sequencing (RNA-Seq) was downloaded from TCGA. The dataset included the expression levels detected for 32,068 different probes on 103 primary melanoma frozen tissues and 368 metastatic melanoma tissues. Among these, 246 DElncRNAs and 856 DEmRNAs were identified (log2fold change > 2, FDR < 0.01) by an analysis in R language. We noted 184 DElncRNAs and 428 DEmRNAs were upregulated in metastatic melanoma, with 62 DElncRNAs and 428 DEmRNAs being downregulated, which were described as a volcano plot (Figures 1A,B).

FIGURE 1

Figure 1. Differentially expressed genes between primary and metastatic melanoma (fold change > 2, FDR < 0.01). (A) Volcano plot of the FDR as a function of weighted fold change for DElncRNAs and (B) DEmRNAs, red dots represent significantly upregulated expressed genes and green dots represent significantly downregulated expressed genes (fold change > 2, FDR < 0.01). (C) Heat map for potential lncRNAs (n = 246) showed significant expression changes, in which 62 were downregulated and 184 were upregulated. Red through green color indicates high to low expression level. (D) Heat map for potential mRNAs showed significant expression changes (n = 856), in which 428 were downregulated and 428 were upregulated.

Hierarchical clustering of the DElncRNA and DEmRNA expression profiles was performed by using “pheatmap” in R language. Heatmaps of the differential expression of the genes allowed for easy discrimination between primary and metastatic melanoma (Figures 1C,D). The top 20 significantly upregulated DElncRNAs are shown in Table 1 and downregulated DElncRNAs in Table 2.

TABLE 1

Table 1. The top 20 upregulated DElncRNAs in metastatic melanoma (p < 0.01).

TABLE 2

Table 2. The top 20 downregulated DElncRNAs in metastatic melanoma (p < 0.01).

The flow diagram of our research design aiming to identify DEGs is shown in Figure 2. We created a gene co-expression network to screen out lncRNAs targeting mRNAs and hub genes among the co-expressed DEGs. Then we performed GO and KEGG analysis on DEGs, followed by analysis of the PPI network. Finally, we verified the prognostic value of DEG expression for patients with melanoma using the comprehensive survival analysis platform Kaplan-Meier plotter. We conducted univariate Cox regression analyses to determine lifetime-predicted values of single genes. We performed ROC curve analyses to obtain the predicted values of multivariate DEGs.

FIGURE 2

Figure 2. Flow diagram of the study design. DElncRNAs and DEmRNAs were investigated between 103 primary and 368 metastatic melanomas (log2fold change > 2, FDR < 0.01)

Co-expression of lncRNAs and mRNAs in Melanoma

A DElncRNA–DEmRNA co-expression network determined by WGCNA was used to reveal the roles and functional mechanisms of lncRNAs and their targeted mRNAs in melanoma (Terai et al., 2016). DElncRNAs and their significantly related mRNAs were imported to draw the network using Cytoscape (version 3.4.3). The co-expression network was composed of 235 network nodes and 762 connection edges, among which 26 upregulated lncRNAs, 9 downregulated lncRNAs, 49 upregulated mRNAs, and 151 downregulated mRNAs were identified to have interregulated relationships (Figure 3). Within the network, several single lncRNAs were co-expressed with multiple mRNAs such as AL51224.1 and AC012313.4, whereas individual mRNAs were targeted by multiple lncRNAs such as AC010266.2, C7orf71, and AC010255.1. Additionally, there were a few co-expression relationships between several lncRNAs or mRNAs, which indicated a highly complex regulatory relationship between lncRNAs and mRNAs in melanoma.

FIGURE 3

Figure 3. lncRNA and mRNA co-expression network in melanoma. The co-expression network established with 235 nodes and 762 edges, in which only edges with weight(w) above a threshold of 0.8 are displayed. The yellow nodes denote upregulated lncRNAs and the blue nodes denote downregulated lncRNAs. The red nodes denote upregulated mRNAs and the green nodes denote downregulated mRNAs.

GO Term Enrichment Analysis

Gene Ontology term enrichment analysis results varied according to GO classification and the expression change of DEmRNAs. GO term enrichment of 363 dysregulated DEGs in the co-expression network was also found among three databases David, GOEAST, GO, and PANTHER. With respect to biological process, the DEmRNAs were significantly enriched in “keratinocyte differentiation,” ”keratinization,” ”peptide cross-linking,” and “epidermis development.” For cellular component, the DEmRNAs were significantly enriched in “extracellular space,” “extracellular region,” “keratin filament,” “desmosome,” “cornified envelope,” “intermediate filament cytoskeleton,” and “intermediate filament.” With regard to molecular function, the DEmRNAs were significantly enriched in “steroid hydroxylase activity,” “structural molecule activity,” ”structural constituent of epidermis,” and “heme binding.” Additional details regarding GO enrichment analysis results are shown in Table 3.

TABLE 3

Table 3. GO term enrichment of dysregulated DEGs in co-expression network coexisting in David, GOEAST, Gene Ontology, and PANTHER.

KEGG Pathway Analysis

Through KEGG pathway analysis, we identified 12 significantly enriched pathways (Table 4): “Steroid hormone biosynthesis” (hsa00140), “Linoleic acid metabolism” (hsa00591), “Salivary secretion” (hsa04970), “Retinol metabolism” (hsa00830), “Aldosterone synthesis and secretion” (hsa04925), “Ovarian steroidogenesis” (hsa04913), “Arachidonic acid metabolism” (hsa00590), “Drug metabolism - cytochrome P450” (hsa00982), “Metabolism of xenobiotics by cytochrome P450” (hsa00980), “Chemical carcinogenesis” (hsa05204), “alpha-Linolenic acid metabolism” (hsa00592), and “Serotonergic synapse” (hsa04726). “Steroid hormone biosynthesis” was also revealed in the GO term analysis. The other results differed from the GO term enrichment analysis, implicating that fairly complicated molecular mechanisms exist in melanoma. Some pathways such as “Hematopoietic cell lineage” and “Ether lipid metabolism” were not significant (P > 0.05).

TABLE 4

Table 4. KEGG pathway enrichment of dysregulated DEGs in melanoma (p < 0.05).

PPI Network Analysis

To predict additional complex relationships, a PPI network of DEmRNAs, consisting of 327 nodes and 456 edges, was constructed using the STRING database and Cytoscape software. The top 18 DEmRNAs with a high degree of connectivity were selected as hub genes for melanoma. The hub genes comprised FGA, FGG, FGFR2, SOX2, S100A9, KRT5, CYP17A1, CYP21A2, FGF1, HSD3B2, FGB, NANOG, FGF10, S100A8, KRT14, FGFBP1, CYP11B2, and CYP11B1. Among these, FGA, FGB and FGG, FGFR2 and FGF1, SOX2 and NANOG, FGFR2 and FGF10, and S100A9 and S100A8 showed the highest node connection degree (Figure 4).

FIGURE 4

Figure 4. The PPI network constructed by the website STRING (https://string-db.org/). A total of 322 proteins are presented with color nodes; interactions are represented with edges. Only edges with weight (w) above a threshold of 0.4 are displayed.

Correlation of DElncRNAs and DEmRNAs With Overall Survival Rate

We next examined whether the DElncRNAs and DEmRNAs were associated with the clinical outcome of patients with melanoma. We performed a univariate Cox regression analysis according to the survival data from TCGA. We performed a univariate Cox regression analysis according to the survival data from TCGA. We found that 14 lncRNAs and 10 mRNAs were significantly associated with patient survival (P < 0.001). Increased expression levels of 12 lncRNAs including AC106882.1, AL365361.1, AC243960.1, AL512631.1, AL928742.1, CHRM3-AS2, LINC00861, LINC01215, LINC01781, LINC02397, LINC02422, and Z98257.1 were correlated with higher survival rate. AC113410.2 and LINC01527 correlated with adverse prognosis. Furthermore, 8 mRNAs SPOCK3, STAP1, THEMIS, TIMD4, TLR10, TMEM156, TNFRSF13B, and TREML2 were associated with favorable outcome in patients with melanoma. Our results suggest that SPRR1B and TENM2, may act as risk factors and independent prognostic indicators. Detailed statistics of Kaplan–Meier survival curves are shown in Figure 5.

FIGURE 5

Figure 5. Kaplan–Meier analysis for overall survival rate of patients from TCGA. Log-rank test was performed to evaluate the survival differences between the two curves (p < 0.001). Overall survival rate for 14 lncRNAs and 10 mRNAs were identified significant.

Cox Regression Analysis and ROC Curve Analysis

To assess the prognostic value of the 246 candidate DElncRNAs and 856 DEmRNAs generated from TCGA database, we conducted Cox regression analysis and time-dependent ROC curve analyses along with calculations of the area under the curve (AUC). In univariate Cox analysis, 486 DEGs were found to have independent significance in predicting clinical outcome (P < 0.05). Among these, 243 DEGs were identified as having great significance (P < 0.001). The top 20 significant predicted-value DEGs are shown in Table 5. NCCRP1, FAM83C, A2ML1, GGT6, TRIM29, RHCG, SPRR2F, PKP1, RHOV, FETUB, KRT17, and CALML3 were identified as independent risk predictors while AL365361.1, AC243960.1, PLA2G2D, TIMD4, PTPRC, DNAJC5B, GPR174, and THEMIS were considered to be positive prognostic factors, which were in accordance with results from the survival rate analysis.

TABLE 5

Table 5. Univariate Cox analysis for lifetime-predicted value of genes (TOP 20).

We further evaluated the potential prognostic value of multivariates among the top 20 DEGs in melanoma using ROC curves. THEMIS, RHCG, PLA2G2D, NCCRP1, AL365361.1, and AC243960.1 were chosen as differences in the multivariate survival curve between high risk and low risk were significant (P = 7.02e-14) (Table 6). The AUC for 3 years was 0.72 (Figure 6A) and that for 5 years was 0.762 (Figure 6B).

TABLE 6

Table 6. Multivariate Cox analysis for lifetime-predicted value of genes coxph(formula = Surv(futime, fustat) ∼ THEMIS + RHCG + PLA2G2D + NCCRP1 + AL365361.1 + AC243960.1, data = rt).

FIGURE 6

Figure 6. The discriminatory ability of the evaluated multivariate DEGs in melanoma was accessed with ROC curve. (A) for 3 year; (B) for 5 years.

Discussion

Metastasis is a crucial factor in the determination of prognosis and survival of melanoma patients. Thus, elucidating the molecular mechanisms underlying melanoma metastasis is urgently needed. lncRNAs are now known to have crucial regulatory functions related to chromatin remodeling and gene expression in the progression of cancers (Wapinski and Chang, 2011). In the current study, a total of 246 DElncRNAs and 856 DEmRNAs were identified. We found 184 DElncRNAs were up-regulated whereas the remainder were down regulated. Furthermore, we investigated the co-expression pattern among 26 up-regulated lncRNAs,9 down-regulated lncRNAs, 49 up-regulated mRNAs and 151 down-regulated mRNAs. “Steroid hormone biosynthesis” and “Linoleic acid metabolism” were the most significant associated pathways according to bioinformatics analysis. In addition, 57 lncRNAs and 36 mRNAs may exhibit significant effects on patient OS rate. Finally, multivariate DEGs in melanoma may have prognostic prediction ability for this disorder.

Emerging evidence indicates that lncRNAs play a vital role in various biological processes such as cell differentiation, immune response, and tumorigenesis (Chen et al., 2018). However, the functions of lncRNAs are not well illustrated (Mercer et al., 2009) and a comprehensive database that provides functions of lncRNAs by experimental researches is still lacking. To further characterize lncRNAs, a “guilt by association” strategy is commonly used to construct a co-expression network of lncRNAs and mRNAs (Liao et al., 2011). In the present study, we used TCGA datasets, which contain evaluated data with high quality and credibility. The transcriptional data was generated by RNA-Seq, which is a revolutionary technology based on next-generation sequencing that is considered the most popular method for studying complete transcriptomes in more detail and with more accurate measurements than microarray and serial analysis of gene expression (SAGE) (Fatica and Bozzoni, 2014). Furthermore, we employed stricter statistical thresholds; fold change > 2 and FDR (q-values) < 0.01 were considered significant to obtain DElncRNAs and DEmRNAs. We systematically investigated DElncRNA-DEmRNA co-expression in a series of statistical models using stringent weight above a threshold of 0.8. This approach identified some new DElncRNAs that showed particularly strong co-expression pairs in our results. For example, AC068985.1 to GABRA6, AC093766.1 and AC007953.2 to GCG, MIR202HG to surrounding mRNAs were predicted to exhibit co-expression relationships (weight > 0.99). In addition, some DElncRNAs, such as AL365361.1, AC243960.1, LINC02397, and LINC02422. AL118505.1, LINC01527, AL512274.1, LINC01215, AC012313.4, FAM41C, SSTR5-AS1, were significant both in co-expression network and predicting survival rate, which may show powerful evidence for the metastasis of melanoma. Notably, the majority of the identified lncRNAs were novel, as only a few lncRNAs have been recorded and characterized in literature to date. To our knowledge, our study is the first to identify novel lncRNAs and mRNAs between primary and metastatic melanoma.

Several previous studies have shown that some lncRNAs involved in the pathogenesis of cutaneous melanoma exert their function through different pathways and that they interact with various molecular targets. Sprouty4-intronic transcript 1 (SPRY4-IT1 or SPRIGHTLY), a cytoplasmically enriched lncRNA derived from an intron of the SPRY4 gene, was the first lncRNA characterized in melanoma in 2011. Khaitan et al. (2011) found that SPRY4-IT1 modulates apoptosis and promotes melanoma proliferation, invasion, and migration in melanoma cell lines. Subsequent studies showed that SPRY4-IT1 can function as an anti-apoptotic factor by binding LPIN2 and can regulate lipid metabolism by avoiding cellular lipotoxicity (Mazar et al., 2014). Patients with melanoma and low levels of SPRY4-IT1 have a longer OS than those with high levels of this lncRNA (Liu et al., 2016). Another study identified the MIR31HG lncRNA gene as showing among the greatest expression changes upon BRAF V600E expression. MIR31HG expression inversely correlated with p16INK4A, and its overexpression may represent a frequent driver event of melanoma progression (Montes et al., 2015). Flockhart et al. (2012) identified BRAF-activated non-coding RNA (BANCR) by analyzing BRAF V600E-mutant human primary melanoma specimens. RNAi-dependent knockdown of BANCR suppressed melanoma cell migration through upregulation of the chemokine CXCL11. It has also been suggested that melanoma proliferation is promoted by BANCR-mediated activation of the ERK1/2 and JNK MAPK pathways (Li et al., 2014). In turn, GAS5 has heretofore been the only lncRNA reported to act as a tumor suppressor in melanoma. Overexpression of GAS5 in the SK-MEL-110 melanoma cell line resulted in a decrease in MMP2 expression involving collagen degradation and a reduction in cell migration (Chen et al., 2016). Furthermore, HOTAIR (Tang et al., 2013), MALAT1 (Tian et al., 2014), UCA1 (Wei et al., 2016), ANRIL (Xie et al., 2016), SAMMSON (Leucci et al., 2016), PVT1 (He et al., 2018), and other lncRNAs are involved in the pathogenesis of malignant melanoma; therefore, they appear to comprise promising biomarkers for early tumor detection or disease progression.

Our GO annotation results indicated that DEmRNAs associated with epidermal cell structure and differentiation might contribute to melanoma tumorigenesis. Kodet et al. (2015) have observed that melanoma cells and neural crest increase expression of keratins 8, 14, 19, and vimentin in co-cultured human primary keratinocytes. FGF-2, CXCL-1, IL-8, and VEGF-A also participate in the activity of melanoma cells on keratinocytes. Expression of keratins 7, 8, 14, 16, 18, and 19 was detected by immunohistochemistry in melanoma. MNF-116, keratins 8, and keratins 18 were higher in metastatic melanoma compared with primary melanoma (Safadi et al., 2016).

The most significant pathway identified by KEGG was steroid hormone biosynthesis. Steroid hormones and their various receptor isoforms, such as androgens and estrogens, play an important role in activating and inhibiting the MAPK and PI3K pathways (de Giorgi et al., 2009, 2011; Mitchell et al., 2014). Notably, the incidence of melanoma increased in men at twice the rate that it has been climbing in women since 1975 (NCI, 2012). It may help explain these phenomena to understand the mechanisms by which estrogen and androgen affect benign and malignant melanocytes. In the genomic pathway, unliganded ERα and ARα in the cell cytoplasm can bind estrogen and androgen, respectively, which leads to dimerization of each receptor. ERα and ARα dimers then traffic to the nucleus and bind to protein co-activators such as CREB to alter proto-oncogene transcription and melanoma proliferation. In the non-genomic pathway, the liganded ERα and ARα activate RAS, which activates BRAF, which activates MEK, which activates ERK, also known as MAPK, which can phosphorylate and therefore activate CREB in the genomic pathway or can directly enhance proto-oncogene transcription and melanoma proliferation (Mitkov et al., 2015).

There are several pathways involved in lipid metabolism including linoleic acid metabolism, arachidonic acid metabolism, and alpha-linolenic acid metabolism. Balogh et al. (2010) observed that cell membrane stress achieved either by heat or benzyl alcohol resulted in pronounced and highly specific alterations in lipid metabolism. A phospholipase C-diacylglycerol lipase-monoacylglycerol lipase pathway was identified in mouse melanoma B16 cells and significantly contributed to the production of several lipid mediators upon stress including arachidonic acid. Schmitt et al. (2015) found that PRDX6 is strongly expressed in most melanoma cells and its expression levels are maintained in a post-transcriptional manner, particularly by EGFR-dependent signaling, while arachidonic acid acts as a key effector of PRDX6-dependent proliferation and inducer of Src family kinase activation. These results subsequently support the biological importance of the emerging evidence of lipid signaling in melanoma.

Another important signaling pathway in our enrichment results is cytochrome P450. Numerous research groups have shown that individuals with polymorphisms that cause defective CYP2D6 are at increased risk for melanoma (Strange et al., 1999). Other genes such as CYP3A4, CYP2C9, ADH4, ADH1A, and CYP1A2 enriched in the cytochrome P450 pathway are observed as metabolically active procarcinogens to genotoxic intermediates, and an association has been found between CYP enzyme activity and the risk to develop several forms of cancer. Amann et al. (2015) found that the retinoic acid receptor (RAR) regulated enzyme CYP26a1 showed a significantly lower expression in the lecithin retinol acyltransferase (LRAT)-overexpressing murine melanoma B16F10 cell line. The absence of LRAT-catalyzed retinol esterification is important for mediating retinoid sensitivity in murine melanoma cells. The potentially complicated mechanisms of the association between cytochrome P450 and retinol metabolism deserve further studies.

It is noteworthy that a single lncRNA was co-expressed with multiple mRNAs and that multiple lncRNAs could be co-expressed with individual mRNAs. Recent studies have shown that most lncRNAs in humans are produced from divergently transcribed protein-coding genes and that the divergent lncRNA/mRNA pairs exhibit coordinated varies in transcription (Sigova et al., 2013; Cabili et al., 2015), representing cis-acting co-expression of lncRNAs and neighboring mRNAs. LncRNAs may also involve in gene expression by targeting distant (trans-acting) coding genes (Cesana et al., 2011; Chu et al., 2011). We identified both neighboring and distant-acting lncRNA and mRNA co-expression. These phenomena suggest highly complex regulatory relationships between lncRNAs and mRNAs in differential co-expression networks.

This study still has limitations. First, the study power is limited by the lopsided size of the samples. Further researches with lopsided larger sample sizes are warranted to validate the reported findings. However, we took the use of multiple kinds of tests into consideration and applied a strict significance threshold to prevent false positive results. Also, we integrated the expression levels of mRNAs and lncRNAs into co-expression profiles to study their characteristics in primary and metastatic melanoma. Experimental validation of the individual roles of lncRNAs is needed. Furthermore, we classified melanoma into primary and metastatic groups because of limited information, although it is likely that melanoma-associated lncRNAs and lncRNA-mRNA co-expression regulation are specific to more detailed clinical phases as well as histological and molecular subtypes. It is essential to perform subgroup analyses by melanoma subtypes in future studies. Finally, the associations of gene expression levels between lncRNAs and mRNAs such as analyzed in this study may not represent the only possible mechanism by which lncRNAs may function in gene regulation; other mechanisms such as promoter demethylation, chromatin remodeling, microRNA silencing, or serving as molecular scaffolds have been reported (Cesana et al., 2011; Chu et al., 2011).

Conclusion

Our work reported the novel dysregulated lncRNAs involved in melanoma determined by using bioinformatics analysis, basing on the large scale of clinical data about lncRNAs between primary and metastatic melanoma. These results might serve as the foundation for establishing future diagnostic or prognostic biomarkers and illuminating potential molecular mechanisms in melanoma. Future studies should be conducted to confirm our findings and to elucidate the detailed tumorigenesis mechanisms of the specific identified lncRNAs.

Author Contributions

LY designed the study, carried out analyses, generated figures, and wrote the manuscript. ZG and SW contributed to statistical assessment, writing, and review. LS designed the study, contributed to writing and review. RT and PL contributed to writing and review.

Funding

This research was funded by the Foundation for Scientific Initiative Plan of Southern Medical University (QD2018N021).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

1. ^ https://cancergenome.nih.gov/
2. ^ https://string-db.org/
3. ^ https://david.ncifcrf.gov/
4. ^ http://omicslab.genetics.ac.cn/GOEAST/
5. ^ http://www.pantherdb.org/
6. ^ http://www.genome.ad.jp/kegg/
7. ^ https://string-db.org/

References

Amann, P. M., Czaja, K., Bazhin, A. V., Rühl, R., Eichmüller, S. B., Merk, H. F., et al. (2015). LRAT overexpression diminishes intracellular levels of biologically active retinoids and reduces retinoid antitumor efficacy in the murine melanoma B16F10 cell line. Skin Pharmacol. Physiol. 28, 205–212. doi: 10.1159/000368806

PubMed Abstract | CrossRef Full Text | Google Scholar

Balogh, G., Péter, M., Liebisch, G., Horváth, I., Török, Z., Nagy, E., et al. (2010). Lipidomics reveals membrane lipid remodelling and release of potential lipid mediators during early stress responses in a murine melanoma cell line. Biochim. Biophys. Acta 1801, 1036–1047. doi: 10.1016/j.bbalip.2010.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabili, M. N., Dunagin, M. C., McClanahan, P. D., Biaesch, A., Padovan-Merhar, O., Regev, A., et al. (2015). Localization and abundance analysis of human lncRNAs at single-cell and single-molecule resolution. Genome Biol. 16:20. doi: 10.1186/s13059-015-0586-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabili, M. N., Trapnell, C., Goff, L., Koziol, M., Tazon-Vega, B., Regev, A., et al. (2011). Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 25, 1915–1927. doi: 10.1101/gad.17446611

PubMed Abstract | CrossRef Full Text | Google Scholar

Cathcart, P., Lucchesi, W., Ottaviani, S., De Giorgio, A., Krell, J., Stebbing, J., et al. (2015). Non-coding RNAs and the control of signalling via nuclear receptor regulation in health and disease. Best Pract. Res. Clin. Endocrinol. Metab. 29, 529–543. doi: 10.1016/j.beem.2015.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Cesana, M., Cacchiarelli, D., Legnini, I., Santini, T., Sthandier, O., Chinappi, M., et al. (2011). A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell 2011, 358–369. doi: 10.1016/j.cell.2011.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., Yang, H., Xiao, Y., Tang, X., Li, Y., Han, Q., et al. (2016). Lentiviral-mediated overexpression of long non-coding RNA GAS5 reduces invasion by mediating MMP2 expression and activity in human melanoma cells. Int. J. Oncol. 48, 1509–1518. doi: 10.3892/ijo.2016.3377

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Sun, Y. Z., Guan, N. N., Qu, J., Huang, Z. A., Zhu, Z. X., et al. (2018). Computational models for lncRNA function prediction and functional similarity calculation. Brief. Funct. Genomics 18, 58–82. doi: 10.1093/bfgp/ely031

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Yan, C. C., Zhang, X., and You, Z. H. (2017). Long non-coding RNAs and complex diseases: from experimental results to computational models. Brief. Bioinform. 18, 558–576. doi: 10.1093/bib/bbw060

PubMed Abstract | CrossRef Full Text | Google Scholar

Chu, C., Qu, K., Zhong, F. L., Artandi, S. E., and Chang, H. Y. (2011). Genomic maps of long noncoding RNA occupancy reveal principles of RNA-chromatin interactions. Mol. Cell 44, 667–678. doi: 10.1016/j.molcel.2011.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

de Giorgi, V., Gori, A., Alfaioli, B., Papi, F., Grazzini, M., Rossari, S., et al. (2011). Influence of sex hormones on melanoma. J. Clin. Oncol. 29, e94–e95. doi: 10.1200/JCO.2010.33.1876

PubMed Abstract | CrossRef Full Text | Google Scholar

de Giorgi, V., Mavilia, C., Massi, D., Gozzini, A., Aragona, P., and Tanini, A. (2009). Estrogen receptor expression in cutaneous melanoma: a real-time reverse transcriptase-polymerase chain reaction and immunohistochemical study. Arch. Dermatol. 145, 30–36. doi: 10.1001/archdermatol.2008.537

PubMed Abstract | CrossRef Full Text | Google Scholar

Fatica, A., and Bozzoni, I. (2014). Long non-coding RNAs: new players in cell differentiation and development. Nat. Rev. Genet. 15, 7–21. doi: 10.1038/nrg3606

PubMed Abstract | CrossRef Full Text | Google Scholar

Flockhart, R. J., Webster, D. E., Qu, K., Mascarenhas, N., Kovalski, J., Kretz, M., et al. (2012). BRAFV600E remodels the melanocyte transcriptome and induces BANCR to regulate melanoma cell migration. Genome Res. 22, 1006–1014. doi: 10.1101/gr.140061.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Gentleman, R., Carey, V., Huber, W., Irizarry, R., and Dudoit, S. (2006). Bioinformatics and Computational Biology Solutions Using R and Bioconductor. Berlin: Springer Science & Business Media.

Google Scholar

Ghazalpour, A., Doss, S., Zhang, B., Wang, S., Plaisier, C., Castellanos, R., et al. (2006). Integrating genetic and network analysis to characterize genes related to mouse weight. PLoS Genet. 2:e130. doi: 10.1371/journal.pgen.0020130

PubMed Abstract | CrossRef Full Text | Google Scholar

He, R. Q., Qin, M. J., Lin, P., Luo, Y. H., Ma, J., Yang, H., et al. (2018). Prognostic significance of LncRNA PVT1 and its potential target gene network in human cancers: a comprehensive inquiry based upon 21 cancer types and 9972 cases. Cell Physiol. Biochem. 46, 591–608. doi: 10.1159/000488627

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang da, W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaitan, D., Dinger, M. E., Mazar, J., Crawford, J., Smith, M. A., Mattick, J. S., et al. (2011). The melanoma-upregulated long noncoding RNA SPRY4-IT1 modulates apoptosis and invasion. Cancer Res. 71, 3852–3862. doi: 10.1158/0008-5472.CAN-10-4460

PubMed Abstract | CrossRef Full Text | Google Scholar

Kodet, O., Lacina, L., Krejèí, E., Dvoøánková, B., Grim, M., Štork, J., et al. (2015). Melanoma cells influence the differentiation pattern of human epidermal keratinocytes. Mol. Cancer 14:1. doi: 10.1186/1476-4598-14-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Leucci, E., Vendramin, R., Spinazzi, M., Laurette, P., Fiers, M., Wouters, J., et al. (2016). Melanoma addiction to the long non-coding RNA SAMMSON. Nature 531, 518–522. doi: 10.1038/nature17161

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Zhang, L., Jia, L., Duan, Y., Li, Y., Bao, L., et al. (2014). Long non-coding RNA BANCR promotes proliferation in malignant melanoma by regulating MAPK pathway activation. PLoS One 9:e100893. doi: 10.1371/journal.pone.0118728

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Q., Liu, C., Yuan, X., Kang, S., Miao, R., Xiao, H., et al. (2011). Large-scale prediction of long non-coding RNA functions in a coding-non-coding gene co-expression network. Nucleic Acids Res. 39, 3864–3878. doi: 10.1093/nar/gkq1348

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Shen, S. K., Xiong, J. G., Xu, Y., Zhang, H. Q., Liu, H. J., et al. (2016). Clinical significance of long noncoding RNA SPRY4-IT1 in melanoma patients. FEBS Open Bio 6, 147–154. doi: 10.1002/2211-5463.12030

PubMed Abstract | CrossRef Full Text | Google Scholar

Mazar, J., Zhao, W., Khalil, A. M., Lee, B., Shelley, J., Govindarajan, S. S., et al. (2014). The functional characterization of long noncoding RNA SPRY4-IT1 in human melanoma cells. Oncotarget 5, 8959–8969. doi: 10.18632/oncotarget.1863

PubMed Abstract | CrossRef Full Text | Google Scholar

Mercer, T. R., Dinger, M. E., and Mattick, J. S. (2009). Long non-coding RNAs insights into functions. Nat. Rev. Genet. 10, 155–159. doi: 10.1038/nrg2521

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitchell, D. L., Fernandez, A. A., Garcia, R., Paniker, L., Lin, K., Hanninen, A., et al. (2014). Acute exposure to ultraviolet-B radiation modulates sex steroid hormones and receptor expression in the skin and may contribute to the sex bias of melanoma in a fish model. Pigment Cell Melanoma Res. 27, 408–417. doi: 10.1111/pcmr.12213

PubMed Abstract | CrossRef Full Text | Google Scholar

Mitkov, M., Joseph, R., and Copland, J. III (2015). Steroid hormone influence on melanomagenesis. Mol. Cell. Endocrinol. 417, 94–102. doi: 10.1016/j.mce.2015.09.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Montes, M., Nielsen, M. M., Maglieri, G., Jacobsen, A., Højfeldt, J., Agrawal-Singh, S., et al. (2015). The lncRNA MIR31HG regulates p16(INK4A) expression to modulate senescence. Nat. Commun. 6:6967. doi: 10.1038/ncomms7967

PubMed Abstract | CrossRef Full Text | Google Scholar

NCI (2012). Surveillance, Epidemiology, and End Results Program (SEER). Washington, DC: National Academic Press.

Google Scholar

Reddy, S. M., Reuben, A., and Wargo, J. A. (2016). Influences of BRAF inhibitors on the immune microenvironment and the rationale for combined molecular and immune targeted therapy. Curr. Oncol. Rep. 18:42. doi: 10.1007/s11912-016-0531-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Safadi, R. A., Bader, D. H., Abdullah, N. I., and Sughayer, M. A. (2016). Immunohistochemical expression of keratins 6, 7, 8, 14, 16, 18, 19, and MNF-116 pancytokeratin in primary and metastatic melanoma of the head and neck. Oral Surg. Oral Med. Oral Pathol. Oral Radiol. 121, 510–519. doi: 10.1016/j.oooo.2015.11.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitt, A., Schmitz, W., Hufnagel, A., Schartl, M., and Meierjohann, S. (2015). Peroxiredoxin 6 triggers melanoma cell growth by increasing arachidonic acid-dependent lipid signalling. Biochem. J. 471, 267–279. doi: 10.1042/BJ20141204

PubMed Abstract | CrossRef Full Text | Google Scholar

Signal, B., Gloss, B. S., and Dinger, M. E. (2016). Computational approaches for functional prediction and characterisation of long noncoding RNAs. Trends Genet. 32, 620–637. doi: 10.1016/j.tig.2016.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Sigova, A. A., Mullen, A. C., Molinie, B., Gupta, S., Orlando, D. A., Guenther, M. G., et al. (2013). Divergent transcription of long noncoding RNA/mRNA gene pairs in embryonic stem cells. Proc. Natl. Acad. Sci. U.S.A. 110, 2876–2881. doi: 10.1073/pnas.1221904110

PubMed Abstract | CrossRef Full Text | Google Scholar

Strange, R. C., Ellison, T., Ichii-Jones, F., Bath, J., Hoban, P., and Lear, J. T. (1999). Cytochrome P450 CYP2D6 genotypes: association with hair colour, Breslow thickness and melanocyte stimulating hormone receptor alleles in patients with malignant melanoma. Pharmacogenetics 9, 269–276. doi: 10.1126/science.1087447

PubMed Abstract | CrossRef Full Text | Google Scholar

Stuart, J. M., Segal, E., Koller, D., and Kim, S. K. (2003). A gene-coexpression network for global discovery of conserved genetic modules. Science 302, 249–255. doi: 10.1126/science.1087447

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, L., Zhang, W., Su, B., and Yu, B. (2013). Long noncoding RNA HOTAIR is associated with motility, invasion, and metastatic potential of metastatic melanoma. Biomed. Res. Int. 2013:251098. doi: 10.1155/2013/251098

PubMed Abstract | CrossRef Full Text | Google Scholar

Terai, G., Iwakiri, J., Kameda, T., Hamada, M., and Asai, K. (2016). Comprehensive prediction of lncRNA-RNA interactions in human transcriptome. BMC Genomics 17:12. doi: 10.1186/s12864-015-2307-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, S., Tian, J., Chen, X., Li, L., Liu, Y., Wang, Y., et al. (2017). Identification of subpathway signatures for ovarian cancer prognosis by integrated analyses of high-throughput miRNA and mRNA expression. Cell Physiol. Biochem. 44, 1325–1336. doi: 10.1159/000485492

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, Y., Zhang, X., Hao, Y., Fang, Z., and He, Y. (2014). Potential roles of abnormally expressed long noncoding RNA UCA1 and Malat-1 in metastasis of melanoma. Melanoma Res. 24, 335–341. doi: 10.1097/CMR.0000000000000080

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Wapinski, O., and Chang, H. Y. (2011). Long noncoding RNAs and human disease. Trends Cell Biol. 21, 354–361. doi: 10.1016/j.tcb.2011.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, Y., Sun, Q., Zhao, L., Wu, J., Chen, X., Wang, Y., et al. (2016). LncRNA UCA1-miR-507-FOXM1 axis is involved in cell proliferation, invasion and G0/G1 cell cycle arrest in melanoma. Med. Oncol. 33:88. doi: 10.1007/s12032-016-0804-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H., Rachakonda, P. S., Heidenreich, B., Nagore, E., Sucker, A., Hemminki, K., et al. (2016). Mapping of deletion breakpoints at the CDKN2A locus in melanoma:detection of MTAP-ANRIL fusion transcripts. Oncotarget 7, 16490–16504. doi: 10.18632/oncotarget.7503

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, X., Hu, Z., Feng, Y., Hu, X., Yuan, J., Zhao, S. D., et al. (2015). Comprehensive genomic characterization of long non-coding RNAs across human cancers. Cancer Cell 28, 529–540. doi: 10.1016/j.ccell.2015.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, X., Zheng, H., Tse, G., Chan, M. T., and Wu, W. K. (2018). Long non-coding RNAs in melanoma. Cell Prolif. 26, e12457. doi: 10.1111/cpr.12457

PubMed Abstract | CrossRef Full Text | Google Scholar