Immune Microenvironment Related Competitive Endogenous RNA Network as Powerful Predictors for Melanoma Prognosis Based on WGCNA Analysis

Cutaneous melanoma is the most life-threatening skin malignant tumor due to its increasing metastasis and mortality rate. The abnormal competitive endogenous RNA network promotes the development of tumors and becomes biomarkers for the prognosis of various tumors. At the same time, the tumor immune microenvironment (TIME) is of great significance for tumor outcome and prognosis. From the perspective of TIME and ceRNA network, this study aims to explain the prognostic factors of cutaneous melanoma systematically and find novel and powerful biomarkers for target therapies. We obtained the transcriptome data of cutaneous melanoma from The Cancer Genome Atlas (TCGA) database, 3 survival-related mRNAs co-expression modules and 2 survival-related lncRNAs co-expression modules were identified through weighted gene co-expression network analysis (WCGNA), and 144 prognostic miRNAs were screened out by univariate Cox proportional hazard regression. Cox regression model and Kaplan-Meier survival analysis were employed to identify 4 hub prognostic mRNAs, and the prognostic ceRNA network consisting of 7 lncRNAs, 1 miRNA and 4 mRNAs was established. After analyzing the composition and proportion of total immune cells in cutaneous melanoma microenvironment through CIBERSORT algorithm, it is found through correlation analysis that lncRNA-TUG1 in the ceRNA network was closely related to the TIME. In this study, we first established cutaneous melanoma’s TIME-related ceRNA network by WGCNA. Cutaneous melanoma prognostic markers have been identified from multiple levels, which has important guiding significance for clinical diagnosis, treatment, and further scientific research on cutaneous melanoma.

Introduction

Cutaneous melanoma is the most malignant skin tumor originating from melanocytes with abnormal proliferation and differentiation, and is the leading cause of skin cancer-related deaths (1). The incidence of melanoma is increasing year by year, and it is characterized by high recurrence rate, high mortality rate, and high drug resistance (2, 3). Early melanoma has a higher cure rate after complete resection, and the mortality rate of stage III/IV patients is as high as 70%, and the 5-year survival rate is less than 16% (4). In recent years, targeted therapy and immunotherapy for melanoma have made certain breakthroughs, while low response rates and severe adverse reactions have limited its long-term effects (5, 6). Therefore, seeking in-depth understanding of the pathogenesis and malignant transformation mechanism of melanoma, exploring more biological targets, and making accurate prediction of the prognosis, will bring new hope for overcoming melanoma’s recurrence and resistance, and improving the survival rate of patients.

Tumor immune microenvironment (TIME) plays an important role in tumorigenesis and therapeutic response. Melanoma cells interact with innate and acquired immune systems, resulting in the inactivation of immune cells with tumor monitoring and killing functions, the activation of tumor protective immune cells, and the formation of TIME with abnormal components and functions, which mediates the immune escape and early metastasis of melanoma (7–9). In order to break the abnormal TIME and mobilize autoimmunity to fight against melanoma, immunotherapy represented by PD-1 and CTLA-4 monoclonal antibody have been worked out, some of which partially delayed the progression of tumor. However, there are still more than half of the patients showing rapid resistance within 6-12 months and even dying of treatment-related adverse drug reactions (10). To clarify the complex tumor-immune regulation in the progression and drug resistance of melanoma, and to find new immunetherapy targets is the only way to to discover new treatment for advanced-stage, recurrent and metastatic melanoma.

Competitive endogenous RNA (ceRNA) network is the emerging post-transcriptional regulation theory in recent years: the non-coding microRNAs (miRNAs) in cells can bind to mRNAs with coding function, resulting in inactivation for translation, unstable structure and rapid degradation of binded mRNAs. At the same time, a variety of other non-coding RNAs (ncRNAs), such as circular RNA (circRNA), long non-coding RNA (lncRNA), and pseudogenes,etc. can sponge absorb miRNA, isolating it from mRNA and forming a competing relationship between ncRNAs and mRNAs, and jointly participating in the endogenous RNA regulatory network, affecting gene translation process (11). CeRNA network regulates the development of most tumors including melanoma (12). Many studies have revealed the specific mechanism by which lncRNA sponge adsorbed miRNA to regulate mRNA expression in melanoma, and clarified the role of various lncRNA-related ceRNA networks in melanoma cells’ proliferation, apoptosis, invasion and migration (13, 14). At the same time, benefiting from the development of high-throughput sequencing technology, more and more studies attempt to outline the post-transcriptional regulation map through systematically analyzing melanoma -related ceRNA networks, with the purpose of seeking for tumor-related prognostic markers (15).

Weighted gene co-expression network analysis (WGCNA) is a method of clustering genes according to expression patterns, systematically analyzing the relationship between gene modules and traits, and classifying gene functions (16, 17). WGCNA method has been widely used in the bioinformatics research on prognostic markers and therapeutic targets of tumors and other diseases.

In order to explore the role of TIME related ceRNA network in melanoma, WGCNA method was used in this study to screen the hub mRNA and lncRNA modules that were significantly related to clinical prognosis from the melanoma transcriptome data of The Cancer Genome Atlas (TCGA) database. The cox survival predicting model was established to screen the hub RNAs and construct the prognostic lncRNA-miRNA-mRNA ceRNA network. mRNAs profiles from Gene Expression Omnibus (GEO) database was used as validation set. We further analyzed the infiltrating immune cells in melanoma samples through The Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) and obtained the phenotype and proportion of immune cells in TIME. Finally, the correlation between immune cell components and prognostic ceRNA network was analyzed. Our research results pointed out the possibility of TIME-associated ceRNA serving as prognostic markers for melanoma, and elucidated the mechanism of melanoma occurrence and development in combination with tumor microenvironment and post-transcriptional regulation.

Materials and Methods

RNA Expression Profiles and Clinical Data of Melanoma Patients

We obtained 471 cutaneous melanoma samples’ mRNA and lncRNA expression profiles, and 450 cutaneous melanoma samples’ miRNA expression profiles from the TCGA database (https://portal.gdc.cancer.gov/). 51 melanoma samples’ mRNA expression profileswere obtained from the mRNA microarray GSE98394 (GPL16791 Illumina HiSeq 2500) in GEO database. Meanwhile, the clinical information of melanoma patients was obtained from TCGA. After the collation of clinical data, 352 cutaneous melanoma samples from 350 patients with complete information of survival time, survival status, and TNM stage were selected for mRNAs and lncRNAs WGCNA analysis. MiRNA expression profiles of 440 melanoma samples with both survival time and survival status information were selected for univariate COX analysis.

Construction of WGCNA Co-Expression Modules and Identifying Prognostic RNAs

WGCNA is an analytical method to identify gene co-expression networks based on topological overlap. The clusters (modules) of highly interconnected genes were determined by hierarchical clustering based on the connectivity and covariation coefficients of the genes. Eigengene expression patterns within each module are condensed into a “Module eigengene (ME).” MEs in the same cluster are considered to have high correlation, consistent expression patterns, and similar biological functions, which is helpful to further explore the functions of different clusters (18). We used the “goodSamplesGenes” function in the “WGCNA” R package to check the missing values of gene expression. After excluding mRNAs and lncRNAs with an average expression amount less than 1, genes with top 25% variance were performed WGCNA analysis to construct co-expression modules (19). The power value is an important soft threshold parameter for defining highly positive correlations among genes in the same module. The “WGCNA” R package was employed to test the independence and average connectivity of different modules under different power values, and the power value corresponding to an independence index of R² = 0.8 was selected. The minimum numbers of genes in the mRNA module and lncRNA module were set to 50 and 30 (20). The “WGCNA” R package was used to cluster thousands of mRNAs and lncRNAs into different modules (identified by arbitrary colors). Based on clinical information, the correlation between modules and clinical phenotypes was calculated, and modules related to prognosis characters such as survival status and survival time were screened to obtain survival-related mRNAs and lncRNAs (WGCNA-mRNAs, WGCNA-lncRNAs) with P < 0.05 as the statistical significance threshold.

Preliminary Screening for Prognostic miRNAs

Cox proportional hazard regression models were constructed with survival outcomes and survival time as dependent variables to analyze the impact of various factors on survival (21). With P < 0.05 as the statistical significance threshold, univariate Cox proportional hazard regression was employed to screen for hub miRNAs (unicox-miRNAs) related to overall survival (OS).

mRNAs Function Annotation

Gene ontology (GO) analysis was performed to WGCNA-mRNAs to annotate the biological processes (BP), cellular component (CC), and molecular function (MF) they involved. Meanwhile, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to annotate the signaling pathways associated with these WGCNA-mRNAs. These analyses were performed by the R software package “clusterProfiler”, with P < 0.05 as the threshold of statistical significance. The “ggplot2” R package was used to visualize the GO and KEGG enrichment results (22).

lncRNA-miRNA and miRNA-mRNA Targeting Relationship Prediction

The miRcode database was used to predict the miRNAs (lnc-pre-miRNAs) targeted by WGCNA-lncRNAs (23). The Venn diagram was employed to take intersection of lnc-pre-miRNAs and unicox-miRNAs, and the mRNAs (mi-pre-mRNAs) targeted by common miRNAs were predicted by Targetscan (24) and miRDB (25) database. Intersection of mi-pre-mRNAs and WGCNA-mRNAs was taken.

Hub Prognostic mRNAs Identifying Through Survival Prediction Model

The univariate Cox proportional hazard regression was employed to further verify the common mRNAs between WGCNA-mRNAs and mi-pre-mRNAs, and screen for hub mRNAs whose expression was related to the OS of patients. MRNAs with P < 0.05 were further evaluated by multivariate Cox proportional hazard regression to identify its ability to predict the prognosis of melanoma. The median risk score obtained through multivariate Cox regression was used as the boundary to divide melanoma patients into high- and low-risk groups. Survival analysis of high- and low-risk groups was performed through the “survival” R package and the survival curve was drawn (26). We also obtained all immunohistochemical pictures of the expression of hub genes in normal skin and melanoma samples from the Human Protein Atlas database (HPA) (http://www.proteinatlas.org), using Image J software to measure the gray-scale of all images to quantify the expression levels of these four protein, and analyzed the obtained gray-scale by GraphPad Prism version 8.0.0 for Windows (San Diego, CA, USA, http://www.graphpad.com). Difference analysis between the two groups was performed by Student’s t test and P < 0.05 was considered statistically significant (27, 28).

Construction of Prognostic ceRNA Network

A lncRNA–miRNA–mRNA ceRNA network with hub prognostic mRNAs as the core was constructed based on the targeting relationships. This ceRNA network elucidates the competitive binding of miRNAs by lncRNAs and mRNAs, thereby revealing the complex post-transcriptional regulatory mechanisms in melanoma. Cytoscape v3.7.2 was employed to visualize the ceRNA network (29).

Analysis of TIME Through CIBERSORT

In order to further study the specific mechanism of abnormal ceRNA network to promote tumors and explore the relationship between ceRNA and TIME, CIBERSORT algorithm was used to identify the proportion of 22 immune cells infiltrated in melanoma samples, and the analysis was repeated 100 times to improve the reliability of the results. The immune cell data of the samples with P < 0.05 were selected for the next correlation analysis (30). Pearson correlation analysis was used to verify the correlation between hub RNAs in ceRNA network and immune cell infiltration.

Statistical Analysis

All statistical analyses were performed by R software (v.3.6.1) and the aforementioned packages.

Results

mRNAs in 3 Co-Expression Modules Preliminary Screened as Melanoma Survival-Related mRNAs by WGCNA Construction

The expression profiles of mRNAs in 470 melanoma samples were explored, the clinic information and the therapies data of the patients were shown in Table 1. After preliminary sorting, we constructed the WGCNA co-expression module for 4938 mRNAs from 352 samples. Firstly, mRNAs were hierarchically clustered (Figure 1A), and Power values determining module independence and average connectivity were screened. When the power value is 4, the independence degree >0.8, which can better reflect the scale-free topology of the co-expression network (Figure 1B). Then 7 mRNAs co-expression modules were established (Figure 1C). The number of eigengenes contained in modules ranged from a minimum of 209 (red modules) to a maximum of 1642 (gray modules).

TABLE 1

Table 1 Clinic information and the therapies data of 470 melanoma patients.

FIGURE 1

Figure 1 Screening survival related mRNAs by WGCNA. (A) Hierarchical clustering tree of 352 melanoma sample mRNAs expression patterns. (B) Power value screening by WGCNA. When the power value is 4, R²>0.8, the average connection degree <100. (C) Clustering and merging of mRNAs co-expression modules. (D) Correlation heatmap of module genes and clinical traits. Red represents positive correlation and blue represents negative correlation. The correlation increases as the the color darkens; (E) Scatterplot of correlation between modules and clinical traits. Brown, blue, and red mRNAs modules are positively correlated with survival status and survival time.

According to the clinical information from TCGA database, the correlation between the module and clinical characters was analyzed (Figure 1D). It was found that brown module, blue module, and red module were highly correlated with survival time and survival status (P < 0.05, Figure 1E). Meanwhile, the blue module was also highly related to the T stage of melanoma (Supplementary Figure 1A). Therefore, it is preliminarily determined that 1401 mRNAs (WGCNA-mRNAs) in the brown, blue, and red modules are related to melanoma prognosis.

GO and KEGG analysis of 1401 WGCNA-mRNAs showed that they were mainly enriched in biological processes such as translation initiation (GO:0006413), endoplasmic reticulum localization (GO:0072599, GO:0070972), interferon gamma response (GO:0034341), extracellular matrix composition (GO:0043062, GO:0030198) (Figure 2A and Table 2), and participate in signal pathways such as ribosomes (hsa03010), antigen presentation (hsa04612), phagosomes (hsa04145), cell adhesion molecules (hsa04514), etc., which were closely related to gene translation, immune response and cell positioning (Figure 2B and Table 3). The P values of the above enrichment analysis results were all <0.05.

FIGURE 2

Figure 2 GO (A) and KEGG (B) analyses of WGCNA -mRNAs. Greater enrichment of specific functions is indicated by larger bubbles and longer columns. Red indicates smaller P values and blue indicates larger P values.

TABLE 2

Table 2 Results of GO enrichment analysis of WGCNA-mRNAs.

TABLE 3

Table 3 Results of KEGG enrichment analysis of WGCNA-mRNAs.

249 lncRNAs in Brown Module Were Preliminary Screened Out as Melanoma Survival-Related by WGCNA Construction

The WCGNA co-expression modules of lncRNAs were constructed in the same way as those of mRNAs. After pre-processing the lncRNAs and clustering (Figure 3A), the power value was screened and determined to be 3 (Figure 3B). Blue (288 lncRNAs), brown (146 lncRNAs), gray (71 lncRNAs), turquoise (315 lncRNAs), and yellow (103 lncRNAs) modules were established (Figure 3C). The module–trait correlation analysis revealed that the brown and yellow modules were closely associated with survival time and survival status (P < 0.05, Figures 3D, E). The brown module was also related to the N stage of the melanoma, and the yellow module was related to the T stage (Supplementary Figures 1B, C). We preliminarily determined that 249 lncRNAs (WGCNA-lncRNAs) in these modules were related to melanoma prognosis.

FIGURE 3

Figure 3 Screening survival related lncRNAs by WGCNA. (A) Hierarchical clustering tree of 352 melanoma sample lncRNAs expression patterns. (B) Power value screening by WGCNA. When the power value is 3, R²>0.8, the average connection degree is <100. (C) Clustering and merging of lncRNAs co-expression modules. (D) Correlation heatmap of module genes and clinical traits. Red represents positive correlation and blue represents negative correlation. The correlation increases as the the color darkens. (E) Scatterplot of correlation between modules and clinical traits. Brown and yellow lncRNAs modules are positively correlated with survival status and survival time.

144 prognostic miRNAs Were Preliminary Screened Through Univariate Cox Regression

We selected the miRNAs expression profiles of 440 samples with both survival time and survival states information for univariate proportional hazard regression. After averaging the expression levels of duplicate miRNAs and screening miRNAs with an average expression level greater than 1, we got a total of 936 miRNAs entering Cox regression. According to the statistical criteria of P < 0.05, 144 miRNAs (unicox-miRNAs), which were related to the prognosis of melanoma were screened out (Supplementary Table 1).

Prediction of Competitive Relationships and Preliminary Screening 48 Hub mRNAs

The miRcode database identified 147 lnc-pre-miRNAs targeted by the 249 WGCNA-lncRNAs, the targeting relationship is visualized in Figure 4A. To screen for miRNAs both involved in the broader ceRNA network and associated with melanoma prognosis, we intersected the 147 lnc-pre-miRNAs and 144 unicox-miRNAs to obtain 43 common miRNAs (Figure 4B).

FIGURE 4

Figure 4 Prediction of the target relationship of lncRNAs, miRNA, and mRNA. (A) Target network between WGCNA-lncRNAs and lnc-pre-miRNAs. (B) Venn diagram was employed to take intersection of unicox-miRNAs and lnc-pre-miRNAs, through which 43 common miRNAs were obtained. (C) Target relationship between common miRNAs and mi-pre-mRNAs. (D) Flow chart of target gene prediction.

The Targetscan and miRDB databases were employed to predict the mRNAs targeted by these 43 common miRNAs. A total of 836 mi-pre-mRNAs were obtained (Figures 4C, D). To screen for mRNAs both involved in the ceRNA network and related to prognosis, we further intersected mi-pre-mRNAs with WGCNA-mRNAs to obtain 48 hub mRNAs (Figure 5A).

FIGURE 5

Figure 5 Screening and Identification of Hub Prognostic mRNAs. (A) The intersection of mi-pre-mRNAs and WGCNA-mRNAs yields 48 hub mRNAs. (B) Survival curve of patients in different risk groups. According to cox regression analysis, the four hub prognostic mRNAs were identified, and the patients were divided into high and low risk groups with median risk value as the boundary. Kaplan-Meier survival analysis showed that the high-risk group had low survival rates (P < 0.001); (C) ROC curve was drawn to verify the accuracy of cox survival model, the area under the curve for predicting 5-year survival rate was 0.802, which was with high predictive efficacy. (D) Expression heatmap of FKBP4, RCC1, GBP1, and PRADC1 hub prognostic mRNAs between high and low risk groups. (E) Immunohistochemical images of the four hub prognostic genes in normal skin and melanoma samples from the Human Protein Atlas database, and gray scale analysis results showed that the expressions of FKBP4, RCC1, and PRADC1 were increased in melanoma and expressions of GBP1 was decreased.

Identifying a 4 Hub mRNAs Prognostic Signature by Cox Regression Model

Univariate Cox proportional hazard regression analysis was employed to verify the correlations between expression levels of these 48 common mRNAs and OS of patients, and 13 putative prognostic mRNAs were identified according to the statistical criteria of P < 0.05. These 13 mRNAs were further analyzed by multivariate COX proportional hazard regression analysis, and finally, four hub mRNAs that could be used to predict melanoma prognosis were selected: PRADC1, FKBP4, RCC1, and GBP1 (Table 4). To verify the accuracy of this model, Kaplan-Meier survival analysis of high and low risk groups was conducted and the result showed that the survival rate of the high-risk group was significantly lower than that of the low-risk group (P < 0.001, Figure 5B). We constructed ROC curves based on the risk scores obtained from the multivariate Cox regression. The area under the 5-year survival curve (AUC) was 0.802, indicating high accuracy and specificity (Figure 5C). The heatmap of the four hub mRNAs between the high and low risk groups is shown in Figure 5D.

TABLE 4

Table 4 Identification of hub prognostic mRNAs through univariate cox regression-lasso regression-multivariate cox regression.

We obtained immunohistochemical images of four hub mRNAs translated proteins in normal skin tissues and melanoma tissues from the HPA database. The gray scale quantification of these images showed that the expression of FKBP4 and RCC1 was significantly upregulated in melanoma, while the expression of PRADC1 was slightly upregulated, and the expression of GBP1 was downregulated (P < 0.05, Figure 5E).

Kaplan-Meier survival analysis was conducted for the four hub mRNAs and survival curves were drawn with TCGA and GSE98394 mRNAs data. The results of TCGA data (Figure 6A) showed that patients with high expression of PRADC1 (P = 0.003), RCC1 (P = 0.004), and FKBP4 (P = 0.053) had low survival rates, and patients with low expression of GBP1 had low survival rates (P = 0.007). Meanwhile, the verification of survival analysis in GSE98394 dataset was completely consisting with the TCGA results (Figure 6B). Patients with high expression of PRADC1 (P = 0.002), RCC1 (P = 0.029), and FKBP4 (P = 0.0008) had low survival rates, and patients with low expression of GBP1 had low survival rates (P = 0.0002).

FIGURE 6

Figure 6 Kaplan-Meier survival analysis for the 4 hub mRNAs with TCGA and GSE98394 mRNAs data. (A) Survival curves of four hub mRNAs with TCGA dataset: Patients with high expression of PRADC1 (P = 0.003), RCC1 (P = 0.004), and FKBP4 (P = 0.053) had low survival rate, while patients with low expression of GBP1 had low survival rate (P = 0.007). (B) Survival curves of four hub mRNAs with GSE98394 dataset: Patients with high expression of PRADC1 (P = 0.002), RCC1 (P = 0.029), and FKBP4 (P = 0.0008) had low survival rates, and patients with low expression of GBP1 had low survival rates (P = 0.0002).

The Four Core Prognostic mRNAs Are Highly Correlated With Prognostic Clinical Information

We also analyzed the differences of the 4 hub prognostic mRNAs expression levels between different subtypes with different TNM stages, total melanoma stages, Breslow thickness, ulcerations, neoadjuvant chemotherapy history, interferon treatment history, systemic system treatment history, and radiotherapy history after diagnosis. The results confirmed that the 4 hub prognostic mRNAs differed among subgroups of different disease states and treatment methods. The expression of GBP1 decreased with the increase of tumor depth (P < 0.0001), and the expression of PRADC1 was the highest in thick melanoma subgroup (P = 0.037, Figure 7A). The expression of GBP1 in melanoma samples with ulcers was lower than that in the non-ulcer group (P = 0.0016, Figure 7B). The expression of GBP1 showed a downward trend with the increase of total melanoma stage and T stage (P < 0.01, Figures 7C, D). The expression level of RCC1 showed an upward trend with the increase of N stage (P = 0.014, Figure 7E). The expression level of PRADC1 in samples with M1 stage was higher than samples with M0 stage (P = 0.047, Figure 7F). In the subgroups with different treatment, the expression of RCC1 in samples after radiotherapy was higher than patients who did not received radiotherapy after diagnosis (P = 0.0054, Figure 7G). The expression of FKBP4 was higher in samples previously received neoadjuvant therapy (P = 0.0033, Figure 7H) and systemic therapy (P = 0.0071, Figure 7I). There is no significant difference in the expression of 4 mRNAs between the interferon treatment and non-treatment groups (Figure 7J). These analysis results were consisting with the Cox regression and survival analysis, indicating that the prognostic mRNAs were highly correlation with melanoma outcome and even have the guidance value to treatment choice in some extent.

FIGURE 7

Figure 7 Identification of the difference of hub prognostic mRNAs’ expression between different clinic phenotype and treatment subgroups. (A) The expression of GBP1 (P < 0.0001) and PRADC1 (P = 0.037) differed among different Breslow depth. (B) The expression of GBP1 was downregulated in melanoma samples with ulcers (P = 0.0016). (C, D) The expression of GBP1 showed a downward trend with the increase of total melanoma stage and T stage (both P value <0.01). (E) The expression level of RCC1 showed an upward trend with the increase of N stage (P = 0.014). (F) The expression of PRADC1 was upregulated in M1 stage samples (P = 0.047). (G) The expression of RCC1 was higher in samples after radiotherapy (P = 0.0054). (H, I) The expression of FKBP4 was higher in samples previously received neoadjuvant therapy (P = 0.0033) and systemic therapy (P = 0.0071). (J) There is no significant difference in the expression of four mRNAs between the interferon treatment and non-treatment groups.

Construction of Melanoma Prognostic lncRNA-miRNA-mRNA ceRNA Regulatory Network

We screened out 4 hub prognostic mRNAs: PRADC1, FKBP4, RCC1, and GBP1. According to the predicted miRNA-mRNA targeting relationship, we identified 1 miRNA: hsa-mir-137. According to the predicted lncRNA-miRNA targeting relationship, 7 lncRNAs that can bind to hsa-mir-137 were matched: HCP5, SNHG14, FGD5-AS1, HCG18, MAGI2-AS3, OIP5-AS1, and TUG1 (Figure 8A). The melanoma prognostic ceRNA regulatory network composed of 7 lncRNAs, 1 miRNA, and 5 mRNAs was finally constructed (Figure 8B).

FIGURE 8

Figure 8 Construction of melanoma prognostic ceRNA network. (A) Flowchart of constructing ceRNA network based on four hub prognostic mRNAs. (B) Melanoma prognostic ceRNA regulatory networks for seven lncRNAs, one miRNA, and four mRNAs were constructed.

Prognostic ceRNA Network Was Closely Correlated With TIME Verified by Comprehensive Analysis

CIBERSORT algorithm was used to analyze the infiltrating immune cell components in melanoma samples. In order to fully ensure the accuracy of the results, the analysis was repeated 100 times and only samples with P < 0.05 were adopted for further analysis. The proportion of 22 immune cells in the TIME is shown in Figures 9A, B. Excluding four non-infiltrating cells, Pearson correlation analysis was employed to identify the correlation among the remaining 18 microenvironmental immune cells (Figure 9C). The results proved the M2 macrophages, activated dendritic cells, neutrophils, eosinophilic cells, plasma cells and helper follicular T cells are positively correlated with each other, and M1 macrophages also have a significant positive correlation with resting dendritic cells (Figure 9E and Supplementary Figures 2 and 3). At the same time, the correlation between 12 RNAs in ceRNA network and immune cell population was analyzed, and it was found that lncRNA-TUG1 had a significant positive correlation with resting NK cells (Figures 9D, F). All the above results proved that the ceRNA network constructed in this study had a certain correlation with the TIME.

FIGURE 9

Figure 9 Correlation between melanoma prognosis ceRNA network and TIME. (A) Histogram of infiltration ratio of 22 immune cells in melanoma microenvironment. (B) Heat map of infiltration ratio of 22 immune cells in melanoma microenvironment. (C) Correlation heatmap of infiltrating immune cells. There were significant positive correlations among M2-type macrophages, activated dendritic cells, neutrophils, eosinophilic cells, plasma cells and helper T cells, and M1-type macrophages and resting dendritic cells. (D) Correlation heatmap of ceRNAs and infiltrating immune cells. (E) Pearson analysis verified that M1-type macrophages were positively correlated with resting dendritic cells in TIME. (F) Pearson analysis verified that lnc-TUG1 was positively correlated with NK cells resting in TIME.

The analysis flow and the schematic diagram of this study is shown in Figure 10.

FIGURE 10

Figure 10 Analysis flow chart and the schematic diagram of this study. (A) Analysis flow chart of data obtaining, processing, and verification. (B) Schematic diagram of this study: The hub prognostic mRNAs (PRADC1, FKBP4, RCC1, and GBP1) compete to bind hub miRNA- hsa-mir-137 with hub lncRNAs (HCP5, SNHG14, FGD5-AS1, HCG18, MAGI2-AS3, OIP5-AS1, and TUG1), which formed a ceRNA network to regulate the prognosis of melanoma. At the same time, resting NK cells in the tumor immune microenvironment were closely related to the hub prognostic lncRNA-TUG1, indicating that the post-transcriptional regulation might have the potential to remodel tumor microenvironment.

Discussion

Melanoma is a type of skin malignant tumor which is extremely threatening to human life. Melanoma progresses rapidly and the mortality rate remains high. Although certain success has been reached in immunotherapy and targeted therapy, early removal of local lesions is the effective treatment at present. However, if it is not detected early and surgically removed, melanoma is highly likely to metastasize and has a strong tendency to spread to other parts of the body and causes serious illness and death (31), putting forward higher requirements for early and accurate diagnosis of melanoma based on effective biomarkers (32). At present, the mechanism of melanoma’s pathogenesis, metastasis and recurrence is still unclear, and many researches are devoted to finding new markers as diagnostic and therapeutic targets. The fate of the tumor is determined by both the tumor cells themselves and the microenvironment (33). As a key post-transcriptional regulatory mechanism, ceRNA network has gradually aroused widespread concern for their important roles in melanoma proliferation, apoptosis, and metastasis. Abnormal TIME also promotes immune escape, drug resistance and metastasis of melanoma (34). In this study, we constructed the melanoma prognostic lncRNA-miRNA-mRNA ceRNA network through WGCNA for the first time, and verified the correlation between ceRNA and TIME, excavated The melanoma prognosis-related markers systematically from multi-dimensions of intracell and extra-cell factors and illustrated the pathophysiology of melanoma.

Previous studies have used WGCNA to identify melanoma prognostic mRNAs. Yang et al. obtained the mRNAs expression profiles of melanoma and normal melanocytic nevi from the GEO database, identified the survival-related mRNAs co-expression modules through WGCNA, and screened out 3 hub prognostic mRNAs through survival analysis (STK26, KCNT2, and CASP12) (3). There are also studies constructing melanoma ceRNA networks through identifying the differentially expressed RNAs in melanoma, and analyzing the relationship between immune cell infiltration and melanoma metastasis (15, 35). In any case, our work is the first research to combine WGCNA, Cox regression, ceRNA network, and TIME. First, we use systematically classify gene expression modules and initially screen for survival-related RNAs (Figures 1 to 3). Then, we took the intersection of the pre-RNAs and WGCNA-RNAs to further screen hub RNAs to improve the prediction effectiveness (Figure 4). Finally, Cox proportional hazards regression model is employed to further identify 4 hub mRNAs and establish the hub prognostic ceRNA network (Figures 5 and 6). The above multiple screening processes ensured the accuracy of the hub ceRNA network as a prognostic marker of melanoma. Combining the analysis results of immune cell components in melanoma microenvironment, the prognostic ceRNA network related to melanoma immune microenvironment was outlined (Figures 7 and 8).

The abnormal components and functions of immune cells in TIME mediate the immune escape, drug resistance, and metastasis of tumor cells. Tumor cells and the microenvironment respond to each other, promoting the reorganization of innate and acquired immune cells components in the tumor microenvironment, and forming a protective immune network for the tumor: differentiation of macrophages in TIME is abnormal, the anti-tumor M1-macrophages is reduced and inactivated, and the tumor-promoting M2 type is increased (34). In the TIME, NK cell activated receptor expression decreases, immune function of NK cells declines (36). Dendritic cell formation is reduced and tended to be quiescent (37). Low expression of T cell receptor (TCR) and interleukin-2 receptor (IL2-R) in CD8+ T cells inhibits their immune function, and number of tumor-protective regulatory T cells increases (38). All above plays key roles in the process of tumor progression. This study found that there was a positive correlation between lnc-TUG1 in melanoma and NK cells resting, suggesting that lnc-TUG1 highly expressed in melanoma may mediate the inactivation of immune surveillance and immune clearance.

The four hub mRNAs (FKBP4, RCC1, GBP1, and PRADC1), one miRNA (hsa-mir-137), and seven lncRNAs (HCP5, SNHG14, FGD5-AS1, HCG18, MAGI2-AS3, OIP5-AS1, and TUG1) in ceRNA network play important roles in the development of melanoma and other various tumors and are widely involved in the immune regulation process.

The FKBP4 gene encodes FKBP51 and FKBP52 proteins that can combine with the immunosuppressive protein FK506, which has been confirmed to be involved in the malignant transformation of melanocytes. FKBP51 protein is highly expressed in melanoma and positively correlated with the vertical growth thickness of tumors, which is an independent prognostic indicator of melanoma. Reducing the expression of FKBP51 in melanoma cells can inhibit the cloning formation of tumor cells and reduce the resistance of tumor cells to ionizing radiation (39). What’s more, FKBP4 is widely involved in a variety of immune processes. FKBP51 protein can mediate the inhibition of FK506 on antigen presentation process, and at the same time, its expression is increased in a variety of immune related diseases such as rheumatoid arthritis (40, 41). These reports all support our discovery that the high expression of FKBP4 predicts poor melanoma prognosis, and that FKBP4 is a hub member of TIME-related ceRNA network. RCC1, as a guanine nucleotide exchange factor of Ran GTPase, plays a key role in cell mitosis, nucleocytoplasmic transport, and nuclear membrane assembly. Qiao et al. found that in cervical cancer associated with papillomavirus E7, RCC1 expression was significantly increased, and destroyed the normal G1 cell cycle checkpoint by regulating transcription factor E2F1 to accelerate tumor cell proliferation (42). Cekan et al. also confirmed that RCC1 could accelerate the cell cycle and help cells resist senescence caused by DNA damage (43). Therefore, we believe that the role of RCC1 in melanoma may be related to its important positive regulation function of cell cycle, and its increased expression will promote the disordered proliferation of tumor cells, and even mediate the resistance of tumor cells to treatment. GBP1 encodes ornithine-binding protein, and plays an important role in both immunity and tumor. Studies have shown that melanoma patients with high expression of GBP family, including GBP1, have a good prognosis in the 30-year follow-uperiod, which is consistent with the result in this study that the decline in GBP1 expression predicts poor prognosis of patients (44). Qiu et al. found that inhibiting the expression of GBP1 in macrophages could promote the polarization of macrophages towards a pro-inflammatory phenotype (45). Increased GBP1 expression can inhibit cell proliferation in an inflammatory environment, while down-regulated GBP1 in various tumors such as melanoma, colon cancer, and breast cancer can help protect cells from the damage of inflammatory response (46). All the above studies have confirmed that the hub mRNAs in the prognostic ceRNA network we screened may determine the fate of tumors by regulating melanoma immune microenvironment. However, there is still a lack of research on PRADC1 in tumors or immunity.

In recent years, our understanding of the function of ncRNA has gradually deepened, and a large number of studies have confirmed that miRNA, lncRNA, and other ncRNA play an important role in regulating tumor and immunity. Hsa-mir-137 (miR-137) is a recognized tumor suppressor gene, and a large number of studies have confirmed that the low expression of miR-137 is significantly correlated with the poor prognosis of melanoma. miR-137 can bind to various mRNAs as ceRNA, reduce the translation level of mRNA, and regulate the biological process in which it participates. Mir-137 can inhibit the migration and invasion of melanoma cells 45–47 by targeting various mRNAs such as PIK3R3, TBX3, c-Met, YB1, EZH2, and MITF (47–49), and inhibit the proliferation and promote the apoptosis of melanoma cells by competitive binding to genes such as SLC1A5, GLO1, CDK6, etc. (50–52). At the same time, Lv et al. found that mouse microglia could activate IL-10R1 and inhibit the pro-inflammatory factors TNF-α and IL-6 by up-regulating the expression of miR-137 to reduce the inflammatory response (53). These results suggest that miR-137 also has significant immunoregulatory functions.

As a recognized tumor suppressor gene, lnc-HCP5 can be used as ceRNA to sponge adsorb miR-12 and further promote the expression of RARRES3 gene, thereby inhibiting the growth of melanoma cells. Wei et al. found that the expression of lnc-HCP5 in melanoma decreased and predicted poor prognosis (54). Major histocompatibility complex (MHC) is an important component of the immune response, mainly consisting of MICA and MICB subtypes. The source gene of lnc-HCP5 is located between MICA and MICB genes, which regulates the activities of various immune cells such as B cells, lymphocytes and NK cells, and plays an important role in HIV, AIDS, and various autoimmune diseases (55). Zhao et al. confirmed the role of lnc-SNHG14 in TIME. They found that in diffuse large B lymphoma, increased expression of lnc-SNHG14 could promote ZEB1 expression through sponge adsorption of miR-5590-3p and further activate PD-L1 and inhibit CD8+ T cells to promote the immune escape of tumor cells (56). Chen et al. found that increasing the expression of lnc-FGD5-AS1 could reduce the inflammation of periodontal ligament cells caused by lipopolysaccharide through miR-142-3p/SOCS6/NF-κB ceRNA axis reduction (57). As an immune-related gene, lnc-HCG18 plays an antitumor role in undifferentiated glioma (58). Luan et al. found that the high expression of lnc-OIP5-AS1 in melanoma could be used as an independent prognostic factor to predict the lower survival rate of melanoma patients, and could sponge adsorb miR-217 as ceRNA to promote glutamine decomposition and accelerate the proliferation rate of melanoma cells (59). Increased expression of lnc-TUG1 in melanoma promoted tumor growth and metastasis by regulating the TUG1/miR-129-5p/AEG-1 and TUG1/miR-29c-3p/RGS1 ceRNA regulatory axis (60, 61). At the same time, reducing the expression of lnc-TUG1 can also suppress the production and release of various immune factors through the TUG1/miR-9-5p/NF-κB1/p50 axis, thus reducing the the inflammatory response of multiple sclerosis. Our study also found that lnc-TUG1 was significantly correlated with quiescent NK cells, suggesting that lnc-TUG1 had the potential to become a tumor immunotherapy target. To date, there have been no reports on functions of SNHG14, HCG18, and MAGI2-AS3 in melanoma, but they constitute a precisely regulated ceRNA network in various solid tumors such as colon cancer, bladder cancer, liver cancer, etc., and regulate the malignant behavior of tumor cells (62, 63). The above studies all support that the ceRNA network we constructed regulates the function of immune cells in the microenvironment and the immune response of tumors, and may play a key role in the formation and evolution of TIME. Meanwhile, the TIME-related ceRNA network in this study has also been proved to be involved in the regulation of the occurrence and development of melanoma and has high prognostic value.

In any case, the limitations of our study suggest that there is still a long way to go in using bioinformatics methods to comprehensively reveal the biological processes of tumors. First of all, although the role of most RNAs in ceRNA network in melanoma or other tumors has been verified by multiple studies, we still found that miR-137 and some lncRNAs showed no significant statistical difference in survival analysis, which means that the amount of data obtained from the public database is still limited, and may limit the accuracy of analysis to some extent. Although this study used WGCNA, Cox regression, Kaplan-Meier survival analysis, and multiple screening methods that combine multiple intersections to improve the prediction accuracy of TIME-related ceRNA for melanoma prognosis, it was still limited to a certain extent by the relatively insufficient number of public data sets under the current situation. Secondly, this study only verified the correlation between the proportion of TIME immune cells and ceRNA, and did not discuss the effect of TIME on melanoma prognosis in detail. In the future, large amounts of data from multiple platforms will be needed to further explore the relationship between immune infiltration and heterogeneity of TIME and melanoma. Finally, although the role of most ncRNAs in this study have been reported by several research respectively, the regulatory and competitive relationship between the lncRNA-miRNA-mRNA axis in melanoma is still lack of experimental verification. Besides, we only screened out the relationship between ceRNA and melanoma TIME, but there is still a lack of relevant research on how this RNA regulates immune cells in the microenvironment.

Conclusion

Our research focused on the identification of prognostic genes by WGCNA, combined with cox regression, Kaplan-Meier survival analysis and CIBERSORT cell composition analysis, etc., established a prognostic ceRNA network related to the melanoma TIME for the first time. We revealed the importance of 4 mRNAs (FKBP4, RCC1, GBP 1, and PRADC1), 1 miRNA (hsa-mir-137), and 7 lncRNAs (HCP5, SNHG14, FGD5-AS1, HCG18, MAGI2-AS3, OIP5-AS1, and TUG1) in reshaping the immune microenvironment and regulating the melanoma process, providing new biomarkers for understanding the mechanism of melanoma occurrence and development at multiple levels, and promoting exploration of prognostic and therapeutic targets.

Data Availability Statement

All datasets presented in this study are included in the article/Supplementary Material.

Author Contributions

Conceptualization: YC and CXL. Methodology: CXL, SW, and LJ. Software: XW, QW, and XS. Validation: XW, YRL, and LY. Formal analysis: YC and YS. Investigation: YL, CYL, and CL. Resources: XW and ZW. Writing (original draft preparation): YC and XW. Writing (review and editing): XW and ZW. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by The National Key R&D program of China (2018YFC1106000).

Conflict of Interest

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.577072/full#supplementary-material

References

1. Tracey EH, Vij A. Updates in Melanoma. Dermatol Clin (2019) 37(1):73–82. doi: 10.1016/j.det.2018.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Schadendorf D, van Akkooi ACJ, Berking C, Griewank KG, Gutzmer R, Hauschild A, et al. Melanoma. Lancet (2018) 392(10151):971–84. doi: 10.1016/S0140-6736(18)31559-9

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yang L, Xu Y, Yan Y, Luo P, Chen S, Zheng B, et al. Common Nevus and Skin Cutaneous Melanoma: Prognostic Genes Identified by Gene Co-Expression Network Analysis. Genes (2019) 10(10):747. doi: 10.3390/genes10100747

CrossRef Full Text | Google Scholar

4. Gershenwald JE, Scolyer RA, Hess KR, Sondak VK, Long GV, Ross MI, et al. Melanoma staging: Evidence-based changes in the American Joint Committee on Cancer eighth edition cancer staging manual. CA: Cancer J Clin (2017) 67(6):472–92. doi: 10.3322/caac.21409

PubMed Abstract | CrossRef Full Text | Google Scholar

5. McKean MA, Amaria RN. Multidisciplinary treatment strategies in high-risk resectable melanoma: Role of adjuvant and neoadjuvant therapy. Cancer Treat Rev (2018) 70:144–53. doi: 10.1016/j.ctrv.2018.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Motofei IG. Malignant Melanoma: Autoimmunity and Supracellular Messaging as New Therapeutic Approaches. Curr Treat Options Oncol (2019) 20(6):45. doi: 10.1007/s11864-019-0643-4

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wu T, Dai Y. Tumor microenvironment and therapeutic response. Cancer Lett (2017) 387:61–8. doi: 10.1016/j.canlet.2016.01.043

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Romano G, Kwong LN. miRNAs, Melanoma and Microenvironment: An Intricate Network. Int J Mol Sci (2017) 18(11):2354. doi: 10.3390/ijms18112354

CrossRef Full Text | Google Scholar

9. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat Med (2018) 24(5):541–50. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Larkin J, Chiarion-Sileni V, Gonzalez R, Grob JJ, Cowey CL, Lao CD, et al. Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. N Engl J Med (2015) 373(1):23–34. doi: 10.1056/NEJMoa1504030

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Tay Y, Rinn J, Pandolfi PP. The multilayered complexity of ceRNA crosstalk and competition. Nature (2014) 505(7483):344–52. doi: 10.1038/nature12986

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Qi X, Zhang DH, Wu N, Xiao JH, Wang X, Ma W. ceRNA in cancer: possible functions and clinical implications. J Med Genet (2015) 52(10):710–8. doi: 10.1136/jmedgenet-2015-103334

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xia Y, Zhou Y, Han H, Li P, Wei W, Lin N. lncRNA NEAT1 facilitates melanoma cell proliferation, migration, and invasion via regulating miR-495-3p and E2F3. J Cell Physiol (2019) 234(11):19592–601. doi: 10.1002/jcp.28559

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Schmidt K, Carroll JS, Yee E, Thomas DD, Wert-Lamas L, Neier SC, et al. The lncRNA SLNCR Recruits the Androgen Receptor to EGR1-Bound Genes in Melanoma and Inhibits Expression of Tumor Suppressor p21. Cell Rep (2019) 27(8):2493–507 e4. doi: 10.1016/j.celrep.2019.04.101

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Huang R, Zeng Z, Li G, Song D, Yan P, Yin H, et al. The Construction and Comprehensive Analysis of ceRNA Networks and Tumor-Infiltrating Immune Cells in Bone Metastatic Melanoma. Front Genet (2019) 10:828. doi: 10.3389/fgene.2019.00828

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol (2005) 4:Article17. doi: 10.2202/1544-6115.1128

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Emilsson V, Thorleifsson G, Zhang B, Leonardson AS, Zink F, Zhu J, et al. Genetics of gene expression and its effect on disease. Nature (2008) 452(7186):423–8. doi: 10.1038/nature06758

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kozela E, Juknat A, Gao F, Kaushansky N, Coppola G, Vogel Z. Pathways and gene networks mediating the regulatory effects of cannabidiol, a nonpsychoactive cannabinoid, in autoimmune T cells. J Neuroinflamm (2016) 13(1):136. doi: 10.1186/s12974-016-0603-x

CrossRef Full Text | Google Scholar

19. Ni Y, Zhang Z, Chen G, Long W, Tong L, Zeng J. Integrated analyses identify potential prognostic markers for uveal melanoma. Exp Eye Res (2019) 187:107780. doi: 10.1016/j.exer.2019.107780

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

21. Fisher LD, Lin DY. Time-dependent covariates in the Cox proportional-hazards regression model. Annu Rev Public Health (1999) 20:145–57. doi: 10.1146/annurev.publhealth.20.1.145

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics: J Integr Biol (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

23. Jeggari A, Marks DS, Larsson E. miRcode: a map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics (2012) 28(15):2062–3. doi: 10.1093/bioinformatics/bts344

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife (2015) 4:e05005. doi: 10.7554/eLife.05005

CrossRef Full Text | Google Scholar

25. Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res (2015) 43(Database issue):D146–52. doi: 10.1093/nar/gku1104

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Bland JM, Altman DG. Survival probabilities (the Kaplan-Meier method). BMJ (1998) 317(7172):1572. doi: 10.1136/bmj.317.7172.1572

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Uhlen M, Zhang C, Lee S, Sjostedt E, Fagerberg L, Bidkhori G, et al. A pathology atlas of the human cancer transcriptome. Science (2017) 357(6352):6352. doi: 10.1126/science.aan2507

CrossRef Full Text | Google Scholar

28. Schneider CA, Rasband WS, Eliceiri KW. NIH Image to ImageJ: 25 years of image analysis. Nat Methods (2012) 9(7):671–5. doi: 10.1038/nmeth.2089

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Su G, Morris JH, Demchak B, Bader GD. Biological network exploration with Cytoscape 3. Curr Protoc Bioinf (2014) 47:8.13.1–24. doi: 10.1002/0471250953.bi0813s47

CrossRef Full Text | Google Scholar

30. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Agaimy A, Specht K, Stoehr R, Lorey T, Markl B, Niedobitek G, et al. Metastatic Malignant Melanoma With Complete Loss of Differentiation Markers (Undifferentiated/Dedifferentiated Melanoma): Analysis of 14 Patients Emphasizing Phenotypic Plasticity and the Value of Molecular Testing as Surrogate Diagnostic Marker. Am J Surg Pathol (2016) 40(2):181–91. doi: 10.1097/PAS.0000000000000527

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Balch CM, Gershenwald JE, Soong SJ, Thompson JF, Atkins MB, Byrd DR, et al. Final version of 2009 AJCC melanoma staging and classification. J Clin Oncol (2009) 27(36):6199–206. doi: 10.1200/JCO.2009.23.4799

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Paget S. The distribution of secondary growths in cancer of the breast. 1889. Cancer Metastasis Rev (1989) 8(2):98–101. doi: 10.1016/j.det.2018.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Ostuni R, Kratochvill F, Murray PJ, Natoli G. Macrophages and cancer: from mechanisms to therapeutic implications. Trends Immunol (2015) 36(4):229–39. doi: 10.1016/j.it.2015.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Xu S, Sui J, Yang S, Liu Y, Wang Y, Liang G. Integrative analysis of competing endogenous RNA network focusing on long noncoding RNA associated with progression of cutaneous melanoma. Cancer Med (2018) 7(4):1019–29. doi: 10.1002/cam4.1315

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ashiru O, Boutet P, Fernandez-Messina L, Aguera-Gonzalez S, Skepper JN, Vales-Gomez M, et al. Natural killer cell cytotoxicity is suppressed by exposure to the human NKG2D ligand MICA*008 that is shed by tumor cells in exosomes. Cancer Res (2010) 70(2):481–9. doi: 10.1158/0008-5472.CAN-09-1688

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Gabrilovich DI. Myeloid-Derived Suppressor Cells. Cancer Immunol Res (2017) 5(1):3–8. doi: 10.1158/2326-6066.CIR-16-0297

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Whiteside TL. Immune modulation of T-cell and NK (natural killer) cell activities by TEXs (tumour-derived exosomes). Biochem Soc Trans (2013) 41(1):245–51. doi: 10.1042/BST20120265

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Romano S, D’Angelillo A, Pacelli R, Staibano S, De Luna E, Bisogni R, et al. Role of FK506-binding protein 51 in the control of apoptosis of irradiated melanoma cells. Cell Death Diff (2010) 17(1):145–57. doi: 10.1038/cdd.2009.115

CrossRef Full Text | Google Scholar

40. Matsushita R, Hashimoto A, Tomita T, Yoshitawa H, Tanaka S, Endo H, et al. Enhanced expression of mRNA for FK506-binding protein 5 in bone marrow CD34 positive cells in patients with rheumatoid arthritis. Clin Exp Rheumatol (2010) 28(1):87–90.

PubMed Abstract | Google Scholar

41. Baughman G, Wiederrecht GJ, Campbell NF, Martin MM, Bourgeois S. FKBP51, a novel T-cell-specific immunophilin capable of calcineurin inhibition. Mol Cell Biol (1995) 15(8):4395–402. doi: 10.1128/mcb.15.8.4395

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Qiao L, Zheng J, Tian Y, Zhang Q, Wang X, Chen JJ, et al. Regulator of chromatin condensation 1 abrogates the G1 cell cycle checkpoint via Cdk1 in human papillomavirus E7-expressing epithelium and cervical cancer cells. Cell Death Dis (2018) 9(6):583. doi: 10.1038/s41419-018-0584-z

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Cekan P, Hasegawa K, Pan Y, Tubman E, Odde D, Chen JQ, et al. RCC1-dependent activation of Ran accelerates cell cycle and DNA repair, inhibiting DNA damage-induced cell senescence. Mol Biol Cell (2016) 27(8):1346–57. doi: 10.1091/mbc.E16-01-0025

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Wang Q, Wang X, Liang Q, Wang S, Xiwen L, Pan F, et al. Distinct prognostic value of mRNA expression of guanylate-binding protein genes in skin cutaneous melanoma. Oncol Lett (2018) 15(5):7914–22. doi: 10.3892/ol.2018.8306

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Qiu X, Guo H, Yang J, Ji Y, Wu CS, Chen X. Down-regulation of guanylate binding protein 1 causes mitochondrial dysfunction and cellular senescence in macrophages. Sci Rep (2018) 8(1):1679. doi: 10.1038/s41598-018-19828-7

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Honkala AT, Tailor D, Malhotra SV. Guanylate-Binding Protein 1: An Emerging Target in Inflammation and Cancer. Front Immunol (2019) 10:3139. doi: 10.3389/fimmu.2019.03139

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Qi J, Wang WW, Chen W, Lu WY, Shang AQ. Mechanism of miR-137 regulating migration and invasion of melanoma cells by targeting PIK3R3 gene. J Cell Biochem (2018) 120:8393–400. doi: 10.1002/jcb.28124

CrossRef Full Text | Google Scholar

48. Peres J, Kwesi-Maliepaard EM, Rambow F, Larue L, Prince S. The tumour suppressor, miR-137, inhibits malignant melanoma migration by targetting the TBX3 transcription factor. Cancer Lett (2017) 405:111–9. doi: 10.1016/j.canlet.2017.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Luo C, Tetteh PW, Merz PR, Dickes E, Abukiwan A, Hotz-Wagenblatt A, et al. miR-137 inhibits the invasion of melanoma cells through downregulation of multiple oncogenic target genes. J Invest Dermatol (2013) 133(3):768–75. doi: 10.1038/jid.2012.357

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Hao S, Luo C, Abukiwan A, Wang G, He J, Huang L, et al. miR-137 inhibits proliferation of melanoma cells by targeting PAK2. Exp Dermatol (2015) 24(12):947–52. doi: 10.1111/exd.12812

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Chen X, Wang J, Shen H, Lu J, Li C, Hu DN, et al. Epigenetics, microRNAs, and carcinogenesis: functional role of microRNA-137 in uveal melanoma. Invest Ophthalmol Visual Sci (2011) 52(3):1193–9. doi: 10.1167/iovs.10-5272

CrossRef Full Text | Google Scholar

52. Luo M, Wu L, Zhang K, Wang H, Zhang T, Gutierrez L, et al. miR-137 regulates ferroptosis by targeting glutamine transporter SLC1A5 in melanoma. Cell Death Diff (2018) 25(8):1457–72. doi: 10.1038/s41418-017-0053-8

CrossRef Full Text | Google Scholar

53. Lv ZC, Cao XY, Guo YX, Zhang XD, Ding J, Geng J, et al. MiR-137-5p alleviates inflammation by upregulating IL-10R1 expression in rats with spinal cord injury. Eur Rev Med Pharmacol Sci (2019) 23(11):4551–7. doi: 10.26355/eurrev_201906_18030

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Wei X, Gu X, Ma M, Lou C. Long noncoding RNA HCP5 suppresses skin cutaneous melanoma development by regulating RARRES3 gene expression via sponging miR-12. OncoTarg Ther (2019) 12:6323–35. doi: 10.2147/OTT.S195796

CrossRef Full Text | Google Scholar

55. Kulski JK. Long Noncoding RNA HCP5, a Hybrid HLA Class I Endogenous Retroviral Gene: Structure, Expression, and Disease Associations. Cells (2019) 8(5):480. doi: 10.3390/cells8050480

CrossRef Full Text | Google Scholar

56. Zhao L, Liu Y, Zhang J, Liu Y, Qi Q. LncRNA SNHG14/miR-5590-3p/ZEB1 positive feedback loop promoted diffuse large B cell lymphoma progression and immune evasion through regulating PD-1/PD-L1 checkpoint. Cell Death Dis (2019) 10(10):731. doi: 10.1038/s41419-019-1886-5

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Chen H, Lan Z, Li Q, Li Y. Abnormal expression of long noncoding RNA FGD5-AS1 affects the development of periodontitis through regulating miR-142-3p/SOCS6/NF-kappaB pathway. Artif Cells Nanomed Biotechnol (2019) 47(1):2098–106. doi: 10.1080/21691401.2019.1620256

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wang W, Zhao Z, Yang F, Wang H, Wu F, Liang T, et al. An immune-related lncRNA signature for patients with anaplastic gliomas. J Neuro-oncol (2018) 136(2):263–71. doi: 10.1007/s11060-017-2667-6

CrossRef Full Text | Google Scholar

59. Luan W, Zhang X, Ruan H, Wang J, Bu X. Long noncoding RNA OIP5-AS1 acts as a competing endogenous RNA to promote glutamine catabolism and malignant melanoma growth by sponging miR-217. J Cell Physiol (2019) 234:16609–18. doi: 10.1002/jcp.28335

CrossRef Full Text | Google Scholar

60. Wang Y, Liu G, Ren L, Wang K, Liu A. Long non-coding RNA TUG1 recruits miR29c3p from its target gene RGS1 to promote proliferation and metastasis of melanoma cells. Int J Oncol (2019) 54(4):1317–26. doi: 10.3892/ijo.2019.4699

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Long J, Menggen Q, Wuren Q, Shi Q, Pi X. Long Noncoding RNA Taurine-Upregulated Gene1 (TUG1) Promotes Tumor Growth and Metastasis Through TUG1/Mir-129-5p/Astrocyte-Elevated Gene-1 (AEG-1) Axis in Malignant Melanoma. Med Sci Monitor Int Med J Exp Clin Res (2018) 24:1547–59. doi: 10.12659/msm.906616

CrossRef Full Text | Google Scholar

62. Xu Z, Huang B, Zhang Q, He X, Wei H, Zhang D. NOTCH1 regulates the proliferation and migration of bladder cancer cells by cooperating with long non-coding RNA HCG18 and microRNA-34c-5p. J Cell Biochem (2019) 120(4):6596–604. doi: 10.1002/jcb.27954

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Yin Z, Ma T, Yan J, Shi N, Zhang C, Lu X, et al. LncRNA MAGI2-AS3 inhibits hepatocellular carcinoma cell proliferation and migration by targeting the miR-374b-5p/SMG1 signaling pathway. J Cell Physiol (2019) 234(10):18825–36. doi: 10.1002/jcp.28521

PubMed Abstract | CrossRef Full Text | Google Scholar