Novel Prognostic Markers for Skin Cutaneous Melanoma

Introduction

The global incidence of skin cutaneous melanoma (SKCM) is increasing more rapidly compared to other types of tumors.¹ Approximately 91% of newly diagnosed skin cancer cases are attributed to SKCM, accounting for about 74% of deaths related to skin diseases.² Overall, surgery, chemotherapy, immunotherapy, and radiation are the primary modalities employed in treating SKCM. Effective biomarkers could potentially serve as valuable guides in optimizing SKCM therapy.³ Previous studies have identified specific biomarkers that could serve as prognostic indicators for certain tumors.⁴ GLT8D1,⁵ Cdc42,⁶ CXCR4,⁶ EZH2,⁷ CD20,⁸ MTAP,⁸ COX-2⁸ and PTEN⁸ have been identified as being associated with the prognosis of SKCM. However, studies on these prognostic biomarkers have primarily focused on their tumor-associated functions as reported in research on other tumors, potentially overlooking some valuable biomarkers specific to SKCM.

In recent years, the advancement of high-throughput research and microarray technology has made it more feasible to explore novel biomarkers for the prognosis of SKCM at the molecular level. However, these studies often focus on individual molecules, overlooking the functional networks they may participate in. Weighted Gene Co-expression Network Analysis (WGCNA) offers a solution by constructing gene co-expression networks and identifying hub genes strongly associated with specific traits.^9,10 Therefore, WGCNA is an excellent method for screening prognostic biomarkers in neoplasms.

In this study, we first utilized transcriptomic analysis to identify significantly differentially expressed genes (DEGs) in SKCM. We then conducted WGCNA to pinpoint the genes most associated with SKCM. By intersecting these genes with the differentially expressed ones and conducting survival analysis, we identified five novel prognostic biomarkers (PLAC8, IL4I1, ZNF80, CCR8, CLEC4C) for SKCM. Furthermore, we experimentally validated these findings using melanoma cell lines and confirmed the functionality of the IL4I1 gene through knockdown experiments. Ultimately, our results underscored the protective roles of these novel prognostic biomarkers in SKCM patients.

Methods

RNA Sequencing (RNA-Seq) Data

RNA-seq count data for SKCM samples were obtained from the The Cancer Genome Atlas Program (TCGA) database (https://portal.gdc.cancer.gov/), while RNA-seq count data for normal skin samples were sourced from the GTEx database (https://gtexportal.org/home/datasets). After merging the TCGA and GTEx data, a gene expression matrix was generated and subsequently processed using R.

Ethics Statement

All human sequencing samples and survival analysis data used in this study were obtained from publicly available databases, specifically TCGA and GTEx. All experimental procedures and protocols used in this research underwent thorough review and received approval from the committee of Dermatology Hospital of Southern Medical University (2023173).

DEGs Analysis

The gene expression matrix was normalized using the Limma package (Ritchie et al).¹¹ DEGs were generated by comparing SKCM samples with normal samples. |log₂Fold change|≥2 and p ≤ 0.01 were the parameters to screen DEGs between SKCM samples and normal samples. The heatmap of all samples was generated by the pheatmap package (Kolde et al¹²), and the volcano map was generated by the ggplot2 package (Wickham et al).¹³ ClusterProfiler package (Yu et al¹⁴) was used for gene ontology (GO) analysis, Kyoto encyclopedia of genes and genomes (KEGG) analysis, and gene set enrichment analysis (GSEA). The correlations between DEGs and diseases (Disease Ontology Semantic enrichment (DOSE) analysis were analyzed using the DOSE package (Yu et al).¹⁵

WGCNA Analysis

The gene co-expression networks were constructed using the WGCNA package (Langfelder et al).¹⁶ Soft power was calculated by the WGCNA package.¹⁶ Once soft power was determined, co-expression modules were generated using algorithms and the blockwiseModules function of WGCNA. Subsequently, the expression level of each module, also known as module eigengene (ME), was calculated. Modules with a similarity of over 75% were merged. The correlations between clinical traits and modules were then calculated. Finally, the most relevant genes associated with SKCM were identified within each module.

Identification of Biomarkers

A key module was identified through WGCNA analysis, with the genes within this module considered as hub genes. Subsequently, SKCM data were stratified into pathology stages I, II, III, and IV based on patient pathology stages. Samples without definite pathology stages were excluded from the analysis. Consequently, a total of 412 SKCM samples were categorized into distinct pathology stages, comprising 77 in stage I, 140 in stage II, 171 in stage III, and 24 in stage IV. Differential expression genes (DEGs) between pathology stages I and IV were identified using the aforementioned methods. Then, hub genes and DEGs of stages I versus IV were took intersection to screen out the novel prognostic biomarkers for SKCM. Cytoscape (Shannon et al¹⁷) was employed to show novel biomarkers and partial hub genes. Transcriptional levels of these genes were produced by the ggstatsplot package. Survival analysis of these genes was processed using the GEPIA database (http://gepia.cancer-pku.cn/) (Tang et al,¹⁸) based on the datasets of TCGA. Additionally, univariate and multivariate Cox analyses were performed to assess the contributions of novel biomarkers using the survival and survminer packages. Default statistical methods provided by these packages and databases were employed to determine significance.

Gene Knockdown

The A375 and HEM cell lines were commercially purchased from Wuxi Xinrun Biotechnology. These cell lines were validated through STR sequencing and confirmed to be free of mycoplasma contamination. Transfection of siRNA was performed using Lipofectamine 3000 reagent according to the manufacturer’s instructions. After 24 hours of transfection, cells were lysed with Trizol reagent to extract RNA for reverse transcription, followed by qPCR to validate knockdown efficiency. Cells collected before and after knockdown were stained with Ki-67 and Annexin-V, washed with PBS for 30 minutes, and then analyzed on a flow cytometer to assess proliferation and apoptosis markers.

Statistical Analysis

Bioinformatics analysis was conducted using default statistical methods. The results of quantitative real-time polymerase chain reaction (qPCR) and flow cytometry were analyzed using non-paired t-tests. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Results

Transcriptional Levels Analysis of RNA-Seq Data

A total of 3563 DEGs were identified, consisting of 1649 up-regulated genes and 1914 down-regulated genes. The heatmap illustrating all genes is depicted in Figure 1A, indicating noticeable differences in transcriptional levels between SKCM and normal samples. Figure 2B displays a volcano plot where red dots signify up-regulated genes, and blue dots represent down-regulated genes. It is evident that DEGs comprise only a small fraction of total genes, with down-regulated genes constituting the majority. Subsequently, GO and KEGG analyses were performed to investigate the functions of DEGs. The GO results (Figure 1C) reveal significant enrichment in biological processes such as epidermis development, skin development, regulation of leukocyte activation, and epidermal cell differentiation. Similarly, cellular component analysis identified significant enrichment in the extracellular matrix, plasma membrane protein complex, and cell-cell junction. Moreover, molecular function analysis showed significant enrichment in receptor regulator activity, receptor ligand activity, and channel activity. As illustrated in Figure 1D, KEGG analysis exhibited significant enrichment in pathways including cytokine-cytokine receptor interaction, neuroactive ligand-receptor interaction, cell adhesion molecules (CAMs), and Staphylococcus aureus infection. The results of DOSE analysis (Figure 1E) indicated significant correlations between DEGs and various conditions such as dermatologic disorders, experimental autoimmune encephalomyelitis, atopic dermatitis, and immunologic deficiency syndromes. Furthermore, KEGG-GSEA analysis (Figure 1F) revealed up-regulation of genes related to the chemokine signaling pathway, Epstein-Barr virus infection, and transcriptional misregulation in cancer, while down-regulation was observed in genes associated with the cAMP signaling pathway, estrogen signaling pathway, and tight junctions. These findings unveil the transcriptional differences between SKCM and normal tissues, further analysis of which may identify prognostically relevant biomarkers.

Figure 1 Transcriptional analysis between normal samples and SKCM samples. (A) The heatmap represents the transcriptional levels of all genes. (B) The volcano plot displays all genes, with up-regulated genes labeled in red and down-regulated genes in blue. The x-axis shows log2 (fold change), while the y-axis represents -log10 p-values. (C) Results of GO analysis show enrichment in biological processes, cellular components, and molecular functions. Node size indicates the number of genes in each term, while color intensity represents the adjusted p-value (all terms had an adjusted p-value of less than 0.05). (D) The bubble chart displays results of KEGG analysis, with the top 10 enriched pathways shown based on gene numbers and adjusted p-values (adjust p-value<0.05). Node size corresponds to the number of genes, and color intensity reflects the adjusted p-value. (E) DOSE analysis of DEGs is depicted in the bar graph, displaying the top 20 enriched pathways based on gene numbers and adjusted p-values (adjust p-value<0.05). (F) KEGG-GSEA analysis was conducted to analyze DEGs, with adjusted p-values of enriched pathways all less than 0.05.

Figure 2 WGCNA was performed to analyze DEGs. (A) The cluster dendrogram illustrates the module eigengenes, with initially clustered modules on the left, followed by similar modules clustering on the right. (B) Correlations between modules and SKCM are depicted, showing R-values and p-values. (C) Clustering between module eigengenes and IA is displayed. (D) Protein-protein interaction networks of hub genes are presented, with triangles representing novel biomarkers and circles representing partial hub genes.

WGCNA Analysis and Identification of Novel Biomarkers

Gene co-expression networks were constructed using data from 19,163 genes across all samples. Soft power 18 was selected to construct 30 modules, as illustrated in Figure 2A. After merging similar modules, 13 modules were obtained. Notably, genes in the grey module were not assigned to any specific modules and were thus excluded from the analysis. The correlation analysis between modules and clinical traits, depicted in Figure 2B, highlighted the Lightgreen module, consisting of 512 genes, as the most relevant to SKCM (r=0.6, p<0.05). Clustering results (Figure 2C) further confirmed the strong correlation between the Lightgreen module and SKCM. By intersecting the 512 genes with DEGs identified in pathology stages I versus IV, PLAC8, IL4I1, ZNF80, CCR8, and CLEC4C were identified as novel potential biomarkers for SKCM prognosis. The network diagram in Figure 2D illustrates the connections among hub genes and these novel biomarkers.

Novel Prognostic Biomarkers for SKCM

As depicted in Figure 3A, the transcriptional levels of the novel biomarkers were markedly elevated in SKCM samples compared to normal samples. Particularly, CLEC4C exhibited negligible expression in normal samples, while its expression levels were significantly upregulated in SKCM samples. Figure 3B illustrates the transcriptional levels of the biomarkers across different pathology stages of SKCM (TCGA dataset). All biomarkers displayed higher expression levels in pathology stage I. Although some differences were observed between pathology stages I and IV, these differences did not reach statistical significance. Furthermore, an overall survival analysis was performed to assess the predictive value of the novel biomarkers for SKCM prognosis. The results revealed that higher expression levels of these biomarkers in SKCM patients were associated with relatively favorable prognoses, indicating the protective effects of the novel biomarkers.

Figure 3 Validation of novel biomarkers. (A) Transcriptional levels between normal samples and SKCM are displayed. (B) Transcriptional levels among different pathology stages are depicted. (C) Survival analysis for novel biomarkers was conducted, with the red line representing high expression in SKCM, and the blue line representing low expression in SKCM.

ZNF80 and CCR8 Were Independent Factors for SKCM Prognosis

Next, we investigated whether these genes could serve as independent predictive factors. As demonstrated in Table 1, univariate Cox analysis indicated that elevated expression levels of PLAC8, IL4I1, ZNF80, CCR8, and CLEC4C were associated with a favorable prognosis. Furthermore, multivariate Cox analysis revealed that high expression levels of ZNF80 and CCR8 were independently associated with a favorable prognosis (refer to Table 1 and Figure 4), suggesting that ZNF80 and CCR8 could function as independent prognostic factors for SKCM. In summary, these findings validate the reliability of the novel biomarkers as prognostic indicators.

Table 1 Cox Analyses

Figure 4 Cox analyses. Forest map of multivariate cox analysis. *P <0.05, ***P<0.001.

IL4I1 is Crucial for the Proliferation and Apoptosis of Melanoma Cells

Next, we examined the expression of these marker genes using HEM (human epidermal melanocytes) and A375 (human malignant melanoma) cell lines, both of which are widely used in melanoma research models. The results showed that, except for CCR8, the expression of other genes was significantly lower in the melanoma cell line A375 (Figure 5A). This finding differs from the analysis results of melanoma tissue from patients (Figure 3A), which may be due to the presence of other cell interference in the overall bulk tissue of patients. The downregulation of these protective genes in melanoma cells suggests their potential importance in melanoma. To validate this hypothesis, we first examined the protein level of IL4I1 and found that, consistent with the gene expression level, the IL4I1 protein was significantly reduced in melanoma cells (Figure 5B). Furthermore, we knocked down IL4I1 (Figure 5C). The results demonstrated that the knockdown significantly promoted the proliferation of melanoma cells (Figure 5D) and inhibited their apoptosis (Figure 5E). These results validate the reliability of the protective markers screened in SKCM.

Figure 5 IL4I1 promotes inhibition of melanoma proliferation and induces apoptosis. (A) qPCR validation of gene expression in HEM and A375 cell lines.(B)Protein-level validation of IL4I1. (C) Western blot analysis confirming IL4I1 protein levels after knockdown. (D) Inhibition of melanoma cell proliferation following IL4I1 knockdown, with the rightmost representing the siRNA-treated group and the leftmost the control group. (E) Apoptosis detection in A375 cells before and after IL4I1 knockdown. Each experiment was independently repeated 2–3 times, and statistical analysis was performed using unpaired t-tests on the results above. *P <0.05, **P <0.01, ***P<0.001, ns: Not significant.

Discussion

SKCM patients exhibit 5-year survival rates of 98.3% in pathology stage I, 63.0% in stage II, and 16.6% in stages III and IV, respectively.¹⁹ Numerous studies indicate that treatment for SKCM is more effective in low-risk patients, whereas high-risk SKCM patients benefit significantly from more rigorous follow-up.^20–25 To date, the treatment of cutaneous melanoma has primarily focused on targeting specific genes and pathogenic pathways, namely the BRAF gene, the MEK/MAPK signaling cascade, the NRAS gene, as well as the KIT gene.²⁶ Therefore, accurately predicting risks in SKCM patients is crucial. In this study, we identified five novel biomarkers through bioinformatics analysis that may serve as prognostic indicators for SKCM. Initially, we compared the transcriptional profiles of normal samples with those of SKCM samples to identify differentially expressed genes (DEGs). Subsequent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the DEGs revealed several significantly enriched terms and pathways, suggesting their potential roles in SKCM pathogenesis. Utilizing Weighted Gene Co-Expression Network Analysis (WGCNA), we identified a lightgreen module strongly correlated with SKCM. Through comprehensive analysis, we pinpointed PLAC8, IL4I1, ZNF80, CCR8, and CLEC4C as potential biomarkers for SKCM prognosis.

Previous studies have revealed that PLAC8 is associated with cell apoptosis, cell differentiation, and disease control.^27–31 In clear-cell renal cell carcinoma, overexpression of PLAC8 has been shown to contribute to a poor prognosis and accelerate tumor progression.³² However, the potential mechanisms and functions of PLAC8 in SKCM prognosis remained ambiguous. In our study, we found that overexpression of PLAC8 may lead to a favorable prognosis in SKCM patients. This may be attributed to the potential induction of tumor cell apoptosis by PLAC8, an interpretation that aligns with our observation of reduced PLAC8 expression in advanced stages of SKCM. In cancer patients, tumor cells and tumor-associated immune cells can produce IL4I1.^33,34 Our experimental results reveal that high expression of IL4I1 in melanoma cells has a protective effect for SKCM patients. However, the multiple functions of IL4I1 in SKCM remain to be further studied. We have found that IL4I1 can inhibit melanoma cell proliferation and promote apoptosis, which may represent a pivotal mechanism underlying the protective role of IL4I1 in melanoma. Regarding ZNF80, its underlying functions in SKCM remain unclear. Our findings suggest that ZNF80 may serve as a novel biomarker for SKCM prognoses. Studies on ZNF80 in the context of tumorigenesis are still scarce; however, given its primary involvement in gene transcription regulation, its role in tumors may involve the modulation of gene expression. CCR8, a marker of tissue-resident memory T cells,³⁵ has been reported to be a necessary molecule for melanoma cells entry into lymph node sinus.³⁶ While the functions of CLEC4C in tumors remain unclear, a study has observed higher expression of CLEC4C in benign breast disease compared to breast cancer patients.³⁷ Also, previous studies have found that CLEC4C⁺ cells are more frequently infiltrated in the lesions of mycosis fungoides.³⁸ However, their specific functional role remains unknown. Our finding suggests the potential of CLEC4C as a favorable prognostic marker in SKCM.

Conclusion

In conclusion, this study has successfully identified five novel prognostic biomarkers, offering valuable insights that could potentially guide both the treatment strategies and prognosis assessment for patients with SKCM. These findings may represent a significant step forward in the quest for more effective management of SKCM, enhancing our ability to tailor treatments to individual patients and improve outcomes.