Identification of Genes Potentially Associated with Melanoma Tumorigenesis Through Co-Expression Network Analysis

Introduction

Melanoma represents one of the most malignant skin cancers; although it only accounts for approximately 4% of skin tumors, it causes an overwhelming 75% of skin-cancer deaths worldwide.¹ Epidemiological survey revealed that the incidence of melanoma is increasing year by year, bringing a heavy economic burden to families and society.² At present, while there are numerous treatment strategies, such as surgery, radiotherapy, chemotherapy and targeted therapy,^3–5 aimed at controlling the disease, the death rate from melanoma is still high due to late-stage diagnosis, thus resulting in poor response to treatment strategies.^6,7 Novel melanoma tumorigenesis genes are required to improve the clinical management and increase life expectancy of patients. Therefore, identifying and analyzing putative candidate genes associated with the transition from normal skin to primary melanoma is particularly urgent.

In recent years, advances in microarray technology coupled with bioinformatics analysis have given researchers the ability to detect putative candidate genes associated with several disease occurrence and tumorigenesis.^8–13 By analyzing the gene expression profiles from 2 GEO databases, Wei et al found that RNF128 in melanoma patients was significantly downregulated in both datasets, indicating that RNF128 is a reliable diagnostic and prognostic biomarker for melanoma.¹⁴ Luan et al found that LncRNA, LICN00520, was significantly overexpressed in melanoma tissue. LICN00520 upregulation was negatively correlated with prognosis of melanoma patients. Their in-vitro study demonstrated that LINC00520 overexpression promotes melanoma cells proliferation, invasion and migration, suggesting that LINC00520 may act as a novel melanoma patient survival biomarker or/and therapeutic target.¹⁵ Another research conducted by Su detected those three genes, CXCL8, THBS1 and KIT, may be the key genes for metastasis melanoma and thus may be the novel therapeutic targets for metastasis melanoma in the future.¹⁶ Taken together, data on melanoma tumorigenesis genes can provide important information to guide treatment, pointing the way for future melanoma research.

In this study, three gene expression profiles, GSE15605, GSE46517 and GSE114445, were reanalyzed and utilized to screened or compare DEGs between primary melanoma and normal skin. Further analysis was performed through bioinformatics to identify the overlapping DEGs (ODEGs) among these three databases. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to characterize the important functional annotations and pathways of the ODEGs involved. Also, a PPI network of the ODEGs was constructed. The MCODE and CytoHubba plug-ins were used to screen the hub genes from the ODEGs. Gene Expression Profile Interactive Analysis (GEPIA) was utilized to evaluate the expression and prognosis of the hub genes in tumor samples. Subsequently, we further analyzed the hub genes related to immune cell infiltration in melanoma and explored the targeted drugs that may have therapeutic effects on melanoma. In summary, our results may provide several reliable immune-related candidate genes and drugs that could be employed to target them, which may help to establish a foundation for the treatment of melanoma.

Materials and Methods

Microarray Data

The Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database is a public functional genomics database of array- and sequence-based data.¹⁷ Therefore, based on search term containing melanoma and inclusion criteria: (1) research type was expression profiling by array, (2) organism was restricted to Homo sapiens, (3) samples contained primary melanoma and normal skin tissue, (4) detection of gene level in tissue samples, (6) studies no less than fifteen samples in each profile. Exclusion criteria: (1) organism was from animal sapiens, (2) samples are at the stage of melanoma tissue with distant metastasis, (3) studies contain less than fifteen samples in each profile, (4) clinical data lack survival time and outcome, as well as detection of gene level. Finally, Gene expression data for primary melanoma and normal skin (GSE15605, GSE46517 and GSE114445) were obtained from GEO database. Specifically, there were 46 primary melanoma and 16 normal skin control samples in GSE15605, 31 primary melanoma and 7 normal skin control samples in GSE46517, and 16 primary melanoma and 6 normal skin control samples in GSE114445. In total, we analyzed 93 primary melanoma and 29 normal skin samples in this study. The microarray data of GSE15605, GSE114445, and GSE46517 were based on an Affymetrix U133 PLUS 2.0 array and an Affymetrix U133A microarray chip, respectively.

Identification of DEGs

The DEGs (GSE15605, GSE114445, and GSE46517) between primary melanoma and normal skin were analyzed by GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). GEO2R is an interactive online tool that enables users to identify DEGs in one or more datasets,¹⁸ which compare the gene expression differences between primary melanoma and normal samples with the R language GEOquery package and the limma package. And those with the threshold of |logFC|>1 and P-value <0.05 were identified as significantly differentially expressed. Next, the volcano plot was constructed at http://www.bioinformatics.com.cn, an online platform for data analysis and visualization. The Venn tool was used to obtain the ODEGs in these three databases. The common DEGs in GSE15605, GSE114445, and GSE46517 were filtered by an interactive tool for comparing lists using Venn diagrams (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

GO and KEGG Pathway Enrichment Analysis

GO term enrichment analysis and KEGG pathway analysis among DEG genes were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID 6.8) tool (https://david.ncifcrf.gov/), which is an online gene enrichment tool, providing a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind large list of genes.¹⁹ GO consisted of three items, namely, the enriched biological process (BP), cellular component (CC) and molecular function (MF). P < 0.05 was set as the threshold value. A heatmap was plotted by online tool software http://www.ehbio.com/Cloud_Platform/front/.

Interaction Analysis of DEGs

A PPI network of the DEGs was constructed using STRING (http://string-db.org/) (version 11.0);²⁰ the data were input in the multiple protein mode, with all parameters using default values. The credibility was 0.70. The PPI network was visually analyzed using Cytoscape 3.6.1 software. Cytoscape’s plug-in Molecular Complex Detection (MCODE) (version 1.4.2) is an application that is used to cluster a given network and find densely connected regions.²¹ After the PPI network is drawn, MCODE is used to identify the most important modules in the PPI network. The selection criteria were as follows: MCODE score > 9, node score cutoff = 0.2, degree cutoff = 2, k score = 2, and maximum depth = 100. KEGG and GO enrichment analyses were performed on the genes in this module.

Hub Gene Selection and Analysis

Hub gene selection was performed using cytoHubba, a plug-in application of Cytoscape. CytoHubba computes 11 methods to identify important nodes in PPI networks.²² Genes appearing at least twice in the top 15 results of each computation method were considered hub genes. GO term and KEGG pathway enrichment further analyses were performed on the results. P < 0.05 was identified as the threshold for significant enrichment.

Expression Level of Hub Genes and Survival Analysis

GEPIA (http://gepia.cancer-pku.cn/) is an interactive online server for analyzing the RNA sequencing expression data of 9736 tumors and 8587 normal samples from the TCGA (The Cancer Genome Atlas) and the GTEx (Genotype-Tissue Expression) projects by a standard processing pipeline. We therefore used these databases to assess the hub genes expression and analysis the association between the hub genes and the overall survival (OS), as well as disease-free survival (DFS) of melanoma samples. By GEPIA, box plot graphs and survival curves were constructed to visualize the expression patterns of the hub genes and their association with prognosis. Log rank test was used to analyses the differences between the two groups. P-values <0.05 were considered to be significant.

cBioPortal Database Analysis

The cBioPortal database (www.cbioportal.org) for Cancer Genomics is an open-source resource that allows users to interactively explore multidimensional cancer genomics data sets. Based on 287 SKCM samples in TCGA, the genetic changes of hub genes were acquired via cBioPortal.

TIMER Database Analysis

Tumor Immune Estimation Resource (TIMER) is a comprehensive online resource for systematical analysis of immune infiltrates across diverse cancer types, allowing users to input specific genes to explore the correlation between gene expression and abundance of immune infiltrates (https://cistrome.shinyapps.io/timer/).²³ Therefore, the 9 hub genes were used to analyze particular immune cell subset (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells) infiltrates in SKCM samples through the TIMER database. P-values < 0.05 were considered to be significant.

Transcription Factor Regulation Network

TRRUST (https://www.grnpedia.org/trrust/) is a database of human and mouse transcriptional regulatory networks, contains 8,444 and 6552 transcription factor (TF)-target regulatory relationships of 800 human TFs and 828 mouse TFs, respectively.²⁴ The TF regulation networks of the hub genes were predicted through the TRRUST database.

Drug–Gene Interaction Analysis

The Drug–Gene Interaction Database (DGIdb) (https://www.dgidb.org/) is an online resource that consolidates disparate data sources describing drug-gene interactions and gene druggability, for mining the druggable genome for precision medicine hypothesis generation.²⁵ Next, the 9 hub genes were used to identify existing and/or potential drug substances associated with these genes through the DGIdb database.

Flowchart of Materials and Methods

The flow diagram of the analysis procedure can be found in the supplementary material Figure S1.

Quantitative Real-Time PCR (qPCR)

Total RNA was isolated using Trizol (Vazyme, China) according to the manufacturer’s instructions. Approximately 1 mg of total RNA from each sample was subjected to first-strand cDNA synthesis using a HiScript III 1st Strand cDNA Synthesis Kit (+gDNA wiper) (Vazyme, China). Real-time PCR was conducted using a SYBR Green qPCR Master Mix according to the manufacturer’s instructions (Bimake, USA). The primers for human were as follows: CCL4 forward (F)-CTCCTCATGCTAGTAGCTGCCTTC, reverse (R)-GGTGTAAGAAAAGCAGCAGGCGGT; CCL5 F: 5-GGCAGCCCTCGCTGTCATCCTCA-3, R: 5-CTTGATGTGGGCACGGGGCAGTG-3; CXCL9 F: 5-GAGTGCAAGGAACCCCAGTAG-3, R: 5-GGTGGATAGTCCCTTGGTTGG-3; CXCL10 F: 5-TGGCATTCAAGGAGTACCTCTC-3, R: 5-GGACAAAATTGGCTTGCAGGA-3; GAPDH F: 5-GGAGTCCACTGGCGTCTTCA-3, R: 5-GTCATGAGTCCTTCCACGATACC-3. Relative mRNA expression levels were calculated using the 2^−ΔΔCt method and normalized to GAPDH expression levels.

Tissue Samples

Ten pairs of human melanoma tissues and adjacent normal tissues were all obtained from the Department of Dermatology, Union Hospital, Tongji Medical College (Wuhan, China) during 2020–2021. They were at least 5 cm away from the tumor site in adjacent normal tissues. And they were stored at −80°C for IHC staining. This process had fully informed consent of the patients. And this study was approved by the Institutional Review Board of Huazhong University of Science and Technology.

Immunohistochemistry (IHC)

All tissues were fixed in 4% paraformaldehyde overnight at 4°C, processed, and cut into 5-μm sections. Masked antigens were uncovered by microwaving the sections for 20 min in citrate buffer, pH 6.0. Then, endogenous peroxidases were quenched, and nonspecific antigen-binding was blocked. The slides were incubated with anti-CCL4 antibody (bs-2475R, Bioss, China), anti-CCL5 antibody (bs-1324R, Bioss, China), anti-CXCL9 antibody (bs-2551R, Bioss, China) and anti-CXCL10 (bs-1502R, Bioss, China) at 4°C overnight. After being washed with PBS (pH 7.6), the slides were treated with streptavidin peroxidase-biotinylated immunoglobulin (Thermo, USA), washed again in PBS. Staining was visualized using 3,3′-diaminobenzidine tetrahydrochloride (Servicebio, China). Counterstaining of the tissue sections was performed with hematoxylin.

Results

Identification of ODEGs

The microarray datasets GSE15605, GSE46517 and GSE114445 were obtained from the public database GEO. Overall, 3270 DEGs (including 1446 upregulated and 1824 downregulated), 1093 DEGs (including 513 upregulated and 580 downregulated), and 1241 DEGs (including 645 upregulated and 596 downregulated) were identified in the GSE15605, GSE46517 and GSE114445 datasets, respectively (supplementary material Table S1). Next, the volcano plots of DEGs are shown in Figure 1A–C, and the expression profiles of the top 10 dysregulated genes between primary melanoma and normal skin in each sample were visualized in a heatmap (Figure 1D–F).

Figure 1 Identification and analysis of ODEGs. (A–C) Volcano plot of the DEGs in GSE15605, GSE46517 and GSE114445. Red and green indicate upregulated and downregulated genes (P < 0.05, Log2FC ≥ 1 or ≤ −1, FC, fold-change threshold), respectively. (D–F) Heatmaps of the top 10 dysregulated genes in each sample of the three datasets. The color scale represents the upregulation (red) or downregulation (blue) of expression. (G) Venn diagram of the ODEGs in the three datasets. (H–I) GO enrichment analyses of the upregulated (H) and downregulated (I) ODEGs. (J and K) KEGG pathway enrichment analyses of the upregulated (J) and downregulated (K) ODEGs. Count, number of DEGs.

Based on these DEGs from the three datasets, we detected 295 ODEGs (157 upregulated and 138 downregulated), as illustrated by a Venn diagram (Figure 1G). Furthermore, these ODEGs were analyzed by GO enrichment (Table S1) and demonstrated that the upregulated genes were primarily enriched in the melanin biosynthetic process, the immune response, melanosomes, melanocyte differentiation and the inflammatory response (Figure 1H); the downregulated genes were primarily enriched in desmosomes, epidermis development, extracellular exosomes, keratinocyte differentiation and transcriptional activator activity, and RNA polymerase II core promoter proximal region sequence-specific binding (Figure 1I). Further KEGG pathway (Table S1) analysis of these DEGs primarily revealed enrichment of the upregulated genes in the Toll-like receptor signaling pathway, cytokine-–cytokine receptor interactions, melanogenesis and the chemokine signaling pathway (Figure 1J); the downregulated genes were primarily enriched in signaling pathways regulating the pluripotency of stem cells and pathways in cancer (Figure 1K). These results suggested that melanogenesis-relevant biological behavior, the immune response and several signaling pathways play a significant role in primary melanoma tumorigenesis.

PPI Network Construction and the Hub Genes Analysis

To further explore the interactions among melanoma tumorigenesis related ODEGs, all the ODEGs of three datasets were inputted and analyzed by the online website STRING, with 0.700 (high confidence) as the minimum required interaction score, and the effect of free proteins was removed. Then, a PPI (Protein–Protein Interaction Networks) network was constructed through Cytoscape software. The network consisted of 149 nodes and 510 lines (Figure 2A). Next, the PPI network was further analyzed among the ODEGs through MCODE with a threshold of MCODE score ≥ 8 and nodes ≥ 8; thus, 2 most important modules were extracted from the PPI network, including module 1 (Score = 9, nodes = 9) and module 2 (Score = 8, nodes = 8) (Figure 2B). Next, the most important modules were selected for further analysis and identification of hub genes through CytoHubba, and the top 10 results in each of the 12 methods (ie, MCC, DMNC, MNC, Degree, EPC, Bottleneck, Eccentricity, Closeness, Radiality, Betweenness, Stress, and Clustering Coefficient) in parameters were sorted, and genes appearing at least third in the top 10 results of each computation method and MCODE score ≥ 8 were considered to be hub genes; thus, 9 hub genes were identified, specifically CXCL10, CXCL9, CCL4, CXCL2, CCL5, NPY1R, PTGER3, NMU, CCL27 (Figure 2C). Next, functional analysis of these 9 hub genes indicated that the genes were primarily enriched in cell chemotaxis, the chemokine-mediated signaling pathway, the extracellular region, the extracellular space, chemokine activity and CXCR3 chemokine receptor binding (Figure 2D). KEGG pathway enrichment showed that these genes were primarily involved in the chemokine signaling pathway, cytokine–cytokine receptor interaction, the Toll-like receptor signaling pathway, the cytosolic DNA-sensing pathway and the TNF signaling pathway (Figure 2E). All these results further suggested that the inflammatory response, immune response and transcriptional activity play an important role in primary melanoma tumorigenesis.

Figure 2 PPI networks of ODEGs and hub genes analysis (A) the PPI network of the ODEGs, Red and blue indicate upregulated and downregulated genes, respectively. (B) The two modules extracted from the PPI network with the highest score. (C) The 9 hub genes were analyzed using CytoHubba. (D and E) GO function and KEGG pathway analysis of hub genes.

Verification of the Expression Patterns and Prognostic Values of Hub Genes

Based on the above results, we found that 9 key genes may play an important role in the switching process from normal skin to primary melanoma tumorigenesis. Thus, to further verify them, the expression patterns of 9 genes in melanoma and normal skin samples were evaluated, and the association of the 9 genes with prognosis was determined using the GEPIA database. The results showed that CCL5, CCL4, CXCL10, and CXCL9 were upregulated (Figure 3A) and demonstrated good OS (Figure 3B) in melanoma samples, while NMU, CXCL2, NPY1R, CCL27, and PEGER3 were downregulated (Figure 3A), and only NMU was associated with worse OS (Figure 3B) in melanoma samples. Moreover, CCL5, CCL4, CXCL10, and CXCL9 upregulation in melanoma samples was associated with good DFS (Figure 3C), although no significance of CXCL10 or CXCL9 was observed in DFS. These results were further validated by Oncomine database (supplementary material Figure S2). In addition, we further analyzed the relationship of between hub genes with better OS and clinical stage of melanoma. Our results indicated that CCL4, CCL5, CXCL9 and CXCL10 are closely related to the clinical stage of melanoma (p < 0.05), while NMU is not (p > 0.05). Notably, CCL4, CXCL9 and CXCL10 have a higher expression in stage 0 and I (supplementary material Figure S3). The relationship of these molecules, CCL4, CCL5, CXCL9, and CXCL10, with the BRAF mutation in melanoma was then investigated further using the cBioPortal database, and the results revealed that there was no link between the expression level of CCL4, CXCL9/10 and the BRAF mutation in melanoma, while CCL5 expression was negatively related to BRAF mutation (supplementary material Figure S4). Our data reveal that these chemokines mainly played their roles in primary melanoma.

Figure 3 Validation of the 9 hub genes in TCGA. (A) The expression patterns of 9 hub genes in melanoma vs normal skin samples. (B) OS curves of 9 hub genes in melanoma samples; (C) DFS curves of 9 hub genes in melanoma samples. P-values < 0.05 were considered to be significant.

Immune Infiltration Analysis

Our results revealed that the inflammatory and immune response may play an important role in this process. An increasing body of evidence suggests that the tumorigenesis, progression and metastasis of melanoma are closely associated with tumor microenvironment.^26–30 Thus, the 9 hub genes and immune infiltration levels were further analyzed through the TIMER database to assess the correlation between them in melanoma. The results of this analysis indicated that there was a positive correlation between CXCL2 expression and the infiltration of CD8⁺ T cells (Cor = 0.108, p = 2.44e−02) and neutrophils (Cor = 0.181, p = 1.05e−04; Figure 4B). PTGER3 expression was positively associated with the infiltration of CD4⁺ T cells (Cor = 0.17, p = 3.19e-04), macrophage cells (Cor = 0.233, p = 5.04e-07), neutrophil cells (Cor = 0.12, p = 1.09e-02) and dendritic cells (Cor = 0.1, p = 3.41e-02; Figure 4F), while a correlation between NMU and CCL27 and the immune cell infiltration level was not observed (Figure 4H–I). However, five genes, CCL5, NPY1R, CXCL10, CXCL9, and CCL4, showed a strong positive correlation between their expression and infiltration of the six immune cell subsets (Figure 4A, C–E and G). These results showed that 7 genes are closely associated with the immune infiltration process of melanoma.

Figure 4 Differently expressed prognostic genes are correlated with the level of immune infiltration. The correlation between the differentially expressed prognostic genes and immune cell infiltration of (A–I) CCL5, CXCL2, NPY1R, CXCL10, CXCL9, PTGER3, CCL4, NMU and CCL27.

Transcription Factor Regulation Network

Based on our results, the transcriptional regulation was also enriched during the initiation of melanoma. The dysregulation of numerous signaling pathways can induce changes in mRNA transcription and drive melanoma initiation and progression.³¹ Next, we further identified the potential transcriptional regulatory networks associated with the hub genes through TRRUST (version 2), which is an open database of human transcriptional regulatory networks. We determined that CCL4, CCL5, CXCL10 and CXCL2 were the key targets of the transcription factors RELA and NFKB1, while CCL5 and CXCL10 were the key targets of the transcription factors IRF7, IRF3 and IRF1 (Table 1). RELA is an important member of NF-kB family, which plays a vital role in regulating the transcription of target genes including the notable proinflammatory cytokines, such as CCL4, CCL5, CXCL10 and CXCL2, which promote an antitumor immunity.³²

Table 1 Key TFs of Hub Genes

Drug–Gene Interaction Analysis

Next, further investigation was conducted using 9 hub genes to explore the possible drugs for melanoma therapy through the DGIdb database. Overall, 46 drug–gene interactions were identified, which were associated with 5 genes (CXCL2, NPY1R, CXCL10, PTGER3, and CCL4). These relationships were further visualized through Cytoscape (Figure 5), whereas the genes CCL5, CXCL9, NMU and CCL27 did not show drug–gene interactions.

Figure 5 Drug–hub gene interaction network. Yellow and green ellipses represent hub genes and drugs, respectively.

Validation by qPCR and IHC Staining in vitro

QPCR was performed to validate the mRNA expression of four hub chemokines in melanoma cell lines and normal melanocytes HEM. The results were consistent with Oncomine database, and indicated that they were abundant in the melanoma cell lines with less metastasis (supplementary material Figure S5). In addition, IHC staining showed that elevated CCL4/5 and CXCL9/10 expression in the melanoma tissues than that in the adjacent normal skin tissues, which was consistent with qPCR results (supplementary material Figure S6).

Discussion

In this study, there GEO datasets, GSE15605, GSE46517 and GSE114445, were selected for bioinformatics analysis. GSE15605 dataset identified novel prognostic marker HMGA2, whose overexpression is strongly associated with regional and distant metastases and serves as an independent predictor of disease-free survival and overall survival in melanoma.³³ GSE46517 dataset supported the thesis that integrated genomic and pathological analyses of staged melanomas provide a productive entry point for discovery of melanoma metastases genes through functional validation results and correlative analysis of human tissues.³⁴ And, GSE114445 dataset suggested that the increased immunosuppressive microenvironment may contribute to the proliferation of melanoma cells, and targeting these suppressive immune cell types may represent an alternative avenue for future immunotherapy.³⁵ Gene expression in these datasets was analyzed between primary melanoma and normal skin, and 9 hub genes were detected in all three datasets, and may have a crucial role during melanoma tumorigenesis. Functional annotation indicated that these hub genes were primarily related to cell chemotaxis, the chemokine-mediated signaling pathway and CXCR3 chemokine receptor binding. The enrichment of genes in KEGG pathway showed that they were primarily associated with the chemokine signaling pathway, cytokine–cytokine receptor interaction, the Toll-like receptor signaling pathway and the TNF signaling pathway. This suggests that chemokines and chemokine-mediated signaling pathway may play an important role in tumor progression. Chemokines are a group of small proteins (~8-14 kDa) that induce chemotaxis, regulate cell proliferation and differentiation of various types of leukocytes via interaction with a subset of G protein–coupled receptors,³⁶ whose main function is related to inflammation and immunity.

CXCL10, a ligand of CXC chemokine receptor-3 (CXCR3), can be secreted by multiple cell types, mainly including T helper cells (Th cells), cytotoxic T lymphocytes (CTLs), dendritic cells, macrophages, natural killer cells (NKs) and some cancer cells.³⁷ CXCL10 is frequently expressed in human disorders affecting inflammation, adaptive immunity, angiogenesis and tumorigenesis.³⁸ Compared with controls, treatment with CXCL10 significantly reduced melanoma tumor growth in C57BL/6 mice in vivo and reduced melanoma cell invasion in vitro.³⁹ CXCL10 deficient mice infected with dengue virus showed a higher mortality rate compared to wild-type mice.⁴⁰ CXCL9, another ligand of CXCR3, can inhibit tumor growth by inducing lymphocytic infiltration to the focal sites.⁴¹ CXCL9 deficient cancer cells are more tumorigenic compared with cancer cells expressing CXCL9.⁴¹ Both CXCR3 and CXCL9 deficient mice are susceptible to immune-mediated kidney disease.⁴² All the results suggested that CXCL9/CXCL10/CXCR3 signaling contribute to anticancer cytotoxic activity. Moreover, both CXCR3 ligands were shown to strengthen the response to some immunotherapy drugs, such as pembrolizumab, nivolumab, atezolizumab and ipilimumab, which are against immune checkpoint PD-1 or CTLA-4.⁴³ These data showed the importance of the CXCL9/CXCL10/CXCR3 axis, which created an opportunity to regulate chemokine expression in the process of tumor therapy and prevention. In addition, we found that the downregulated NMU was associated with worse OS, while upregulated CCL4, CCL5, CXCL9 and CXCL10 were related to better OS in melanoma patients, these results indicated that chemokine members (CCL4, CCL5, CXCL9 and CXCL10) and NMU could play important roles in melanoma initiation and development. CXCL9 and CXCL10 are ELR(Glu-Leu-Arg) motif-negative CXC chemokines secreted by multiple cell types, and could ameliorate angiogenesis and tumor growth in various cancers, such as gastric,⁴⁴ prostate,⁴⁵ lung,⁴⁶ cervical cancer⁴⁷ and melanoma,³⁹ suggesting that they are important biomarkers in melanoma tumorigenesis.

Melanoma is frequently infiltrated by inflammatory and immune cells, therefore promoting or inhibiting tumor growth.³⁹ We observed that CXCL2 and PTGER3 was positively associated with partial infiltration of cells, and the genes CCL5, NPY1R, CXCL10, CXCL9 and CCL4 were all positively correlated with the infiltration of 6 immune cell types, including B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils and dendritic cells. In colorectal cancer, CXCL10 was found to be related to CXCR3⁺CD8⁺ T cells and IFN-γ⁺CD4⁺ T cells of the tumor site, potentiating anti-tumor activity through production of Granzyme B and IFN-γ from CD8⁺ T cells and CD4⁺ T cells, respectively.⁴⁸ Treatment CXCL10-Ig to mice engrafted with myeloma cells can significantly inhibit the tumor growth by reducing the number of Treg cells and accumulating CD4⁺, CD8⁺ and NK cells at the tumor site.⁴⁹ Another study showed that CCL5 secreted by cancer cells, following recruiting NK cells to the tumor site, acts its function via paracrine signaling.⁵⁰ CCL4/CCL5 interact with their receptor CCR1/CCR5, which facilitates the recruitment of γδ T cells from peripheral blood and/or peritumoral regions to the tumor sites, thus enhancing the antitumor immunity and improving patients’ prognosis.⁵¹ All the findings suggested that CCL4, CXCL10 and CCL5 may be potent biomarkers or targets for melanoma.

In addition, dysregulation of transcription can drive a wide range of diseases. Our results showed that RELA and NFKB1 could be the key transcription factors of CCL4, CCL5, CXCL10 and CXCL2; IRF7, IRF3 and IRF1 are the key TFs associated with CCL5 and CXCL10. NF-κB is a heterodimer composed of the RELA (p65) and NFKB1 (p50) subunits that belong to the Rel/NF-κB family of transcription factors. Several studies have demonstrated that RELA and NFKB1 play important roles in the control of cell proliferation and apoptosis and, hence, oncogenesis,^52,53 and its activity significantly contributes to the tumorigenicity, angiogenesis, and metastasis of human melanoma cells.⁵⁴ Interferon regulatory factors (IRFs) are a family of homologous proteins that play key roles in mediating signal transduction associated with the immune response and immunomodulation in Homo sapiens.⁵⁵ Accordingly, IRF1 and IRF2 could promote the response of Th1 cells, whereas IRF3 and IRF7 are involved in antibacterial and antiviral innate immunity,⁵⁶ suggesting that they are vital elements of the immune system in disease.

To further understand and find the potent drugs of melanoma, we explored the drug–gene interactions based on the DGIdb database, 46 drugs associated with 5 genes (CXCL2, NPY1R, CXCL10, PTGER3, and CCL4) in melanoma were identified. Oxaliplatin, for instance, is widely used to treat gastrointestinal and other malignancies. Concomitant treatment of anti-PD-1 (an immune-checkpoint inhibitor) with trifluridine/tipiracil (FTD/TPI, a new antimetabolite agent), oxaliplatin can be synergized to promote induced immunogenic cell death (ICD). This enhances the antitumor efficacy of dual chemotherapy in both mice and human colon carcinoma cell lines.⁵⁷ Clinical trials have also revealed that adjuvant oxaliplatin, leucovorin, fluorouracil bolus and fluorouracil infusion improved DFS in patients with rectal cancer after receiving preoperative chemoradiotherapy.⁵⁸ Moreover, beraprost, a chemically stable and orally active prostacyclin analog, is often used to treat chronic arterial occlusive disease. Interestingly, a recent study showed that mice with melanoma were treated with beraprost for 3 weeks, the number and size of lung tumor metastases are significantly reduced. And, beraprost treatment also raised the number of tumor-related pericytes and improved intratumor hypoxia which was validated by scanning electron microscopy and immunohistochemistry.⁵⁹ These data suggest that drug–gene interaction pairs we predicted may be of value, revealing the potential candidate for melanoma treatment.

Conclusion

In our study, we identified 9 hub genes associated with tumorigenesis and progression in melanoma. They were closely associated with several immune cell infiltrates and could be regulated by several transcription factors. Taken together, our findings will provide new understanding for melanoma immunotherapy, assisting the doctors in making decision and selecting efficacious drug. However, further in vitro and in vivo studies should be performed to verify the roles played by these 9 hub genes and analyze their roles in the tumorigenesis of melanoma.