Clinical specimens were provided by Guizhou Medical University Affiliated Hospital (Guiyang, China), comprising 92 pairs of tumor and adjacent samples (1-2 cm from tumor tissue). These samples encompassed 24 trunk subtype samples and 68 other subtype samples. Approval for this study was obtained from the Human Characteristics Ethics Committee of Guizhou Medical University, and the principles of the Helsinki Declaration were strictly followed. Informed consent was also obtained from patients providing samples.

Transcriptomic data and related clinical information were sourced from the TCGA database (https://portal.gdc.cancer.gov/projects/TCGA-SKCM). The RNA-sequencing dataset encompasses 472 samples, including 78 samples from Trunk Cutaneous Melanoma (TCM) and 394 samples from other melanoma subtypes. A threshold of 140, post-normalization, was applied to identify outliers using the hierarchical clustering algorithm. The cut-off for FPKM was established at 0.5 to exclude genes with low expression levels. A threshold of 0.5 was set during gene selection. Following merge and batch normalization, differential expression genes (DEGs) are identified using ‘‘sva’’ R package (13–15). Moreover, DEGs were discerned between tumor and normal tissues using a cut-off value of |Log2 fold-change (FC)| > 0.5 and an adjusted P-value < 0.05 (16).

Genes were analyzed for the construction of a scale-free network via the ‘‘WGCNA’’ methodology (17, 18). This expression network, predicated on weighted gene co-expression data, was established utilizing the top quartile of genes exhibiting the highest variance. Initially, an expression data-derived similarity matrix was generated to compute the absolute values of Pearson correlation coefficients between gene pairs. Subsequently, this similarity matrix was transformed into an adjacency matrix, with the application of a soft threshold that accentuated strong connections while diminishing the significance of weaker correlations within the matrix. This adjacency matrix was further converted into a topological overlap matrix (TOM), designed to more accurately depict the strength of relationships and the degree of connectivity among genes within the dataset. Serving as a measure for assessing network connectivity among genes, TOM was employed as an input for hierarchical clustering analysis. The ‘‘DynamicTreeCut’’ algorithm, utilized within the ‘‘WGCNA’’ R package, facilitated the identification of network modules (19). Following module identification, gene significance (GS) was evaluated, correlating these genes with other biological information. The higher the GS value, the more substantial the prognostic relevance for patients. Hence, an analysis employing Pearson’s correlation coefficient was conducted to ascertain the correlation between identified modules and clinical features, including patient age and cancer stage. Through this approach, ‘‘WGCNA’’ aids in elucidating the connections between gene expression patterns and clinical characteristics, offering vital insights for the discovery of biomarkers and the formulation of therapeutic strategies for disease management.

GO and KEGG enrichment analyses were conducted using the R package ‘‘clusterProfiler’’ to explore pivotal modules, with P < 0.05 defined as significant enrichment (20). These analyses aimed at identifying potential signaling pathways, elucidating the biological processes and pathways crucial in the occurrence and progression of melanoma (21).

PPI network was constructed using the online STRING database (https://string-db.org/). During network construction, disconnected nodes were removed to ensure network connectivity. Subsequently, the network graph of key modules was created using Cytoscape 3.72, revealing crucial interactions by identifying the top 30 ranked genes. Each node in the PPI network represents a protein encoded by a gene, while edges depict interactions between two proteins. To explore the potential of DSC2 and DSC3 as novel targets for regulatory intervention, an analysis was conducted to identify genes that show co-expression with these markers. This analysis highlighted twelve genes exhibiting the highest co-expression coefficients with DSC2 and DSC3. The findings, represented in a circular graph, underscore the importance of further investigation into these associated genes, suggesting they warrant closer scrutiny.

Survival analysis was conducted using the ‘‘survival’’ and ‘‘survminer’’ packages in R language (22). Initially, patients were stratified into high and low expression groups based on the median gene expression levels, and the association between these gene expression levels and prognosis was determined using the Log-rank test.

Subsequently, exploration of subpopulation clustering in SKCM at the single-cell level was carried out utilizing the TISCH2 single-cell database (23). Briefly, the ‘‘Seurat’’ package is employed for scrutinizing scRNA-seq data of SKCM (18). Initially, following the exclusion of cells expressing fewer than 250 or exceeding 6000 genes, a logarithmic normalization is performed on gene expression. Subsequently, employing the functions ‘‘FindNeighbors’’ and ‘‘FindClusters’’, individual cells are clustered into distinct subgroups. To further elucidate the roles of DSC2 and DSC3 in immune infiltration, the ‘‘CIBERSORT’’ algorithm was employed to calculate immune infiltration scores in each sample (24). Samples were dichotomized based on the median gene expression levels, comparing the infiltration of various types of immune cells between these two groups and assessing DSC2 and DSC3 respectively. The correlation analysis between gene expression and immune cell infiltration was visualized using the ‘‘ggplot2’’ package, including box plots and correlation scatter plots. Finally, the correlation between genes and immune checkpoints was analyzed using the ‘‘corrplot’’ package and visualized through heat maps. Furthermore, specific gene associations with immune checkpoint genes were evaluated.

In recent years, immunotherapy has emerged as a promising anti-tumor strategy (25), demonstrating favorable anti-tumor effects in SKCM treatment. Hence, further investigation into the role of oxidative stress-related genes in SKCM immunotherapy was conducted. The association between immunotherapy response and oxidative stress-related genes was explored using TIDE analysis. Additionally, the BEST database illustrated the impact of oxidative stress signaling pathways on chemotherapy drug resistance. Immunohistochemical staining and ROC curve analysis were employed to evaluate the protein expression levels of core genes in SKCM and normal samples, and their diagnostic value was determined through survival analysis (26).

GSVA is a non-parametric unsupervised analytical method used to assess the enrichment of gene sets in transcriptome data (27). The score of oxidative stress were determined by the ‘‘ssGSEA’’ using ‘‘GSVA’’ R package, and the certain signatures that respond to oxidative stress were obtained from database (http://www.informatics.jax.org/vocab/gene_ontology/GO:0006979). This method transforms gene expression data into an expression matrix characterized by a specific set of genes (28). Besides, Pearson correlation analysis was conducted to validate the correlation between hub genes and oxidative stress scores. Through this process, we could assess association between screened genes and oxidative stress.

In an ambient laboratory environment, tissues were subjected to a 4% Paraformaldehyde fixation process. Subsequently, these specimens were encased within paraffin blocks and dissected into sections of 4 µm in thickness. The process of removing paraffin was conducted utilizing xylene at a temperature of 60°C, succeeded by a graded rehydration sequence employing ethanol concentrations of 100%, 80%, 60%, and 40%. To inhibit the intrinsic peroxidase reactions present in the samples, 3% H₂O₂ was utilized, while sodium served the purpose of facilitating antigen retrieval. Following a 16-hour incubation period at 4°C in a medium containing 5% bovine serum albumin, the sections underwent a further incubation with an array of primary antibodies (Santa Cruz, CA, USA), specifically (DSC2; 1:200), (DSC3; 1:200), (DSG1; 1:200), (KRT6B; 1:200), (PKP1; 1:200), and (PKP3; 1:400). Subsequent to this, for a duration of 2 hours at an ambient temperature, the sections were exposed to secondary antibodies targeted against mouse and rabbit immunoglobulins. Visual documentation of the samples was achieved using a high-resolution light microscopy technique, employing magnification levels of x200 and x400. Following this, a staining protocol involving 3,3’ diaminobenzidine and hematoxylin was applied for a minute at ambient temperature. The evaluation of the staining intensity for the genes of interest was performed, with scores assigned as follows: 0 indicating no staining, 1 for weak positivity, 2 for moderate positivity, and 3 for strong positivity, facilitating the quantification of protein expression. The Image Pro Plus software was employed to determine the scoring based on the proportion of positively stained cells, with a scoring rubric of 0-2 for low expression, 3-4 for moderate expression, and 5-6 denoting high expression. P-value < 0.05 was interpreted as indicative of a statistically significant variance.

All statistical analyses were conducted using packages in the R programming language (version 3.6.3). Prior to statistical analysis, normality tests were performed to determine the appropriate statistical methods. Analysis of normal and non-normal data was conducted using unpaired Student’s t-test and Wilcoxon test, respectively (29). Pearson method was employed for correlation analysis, with P < 0.05 considered statistically significant. * P < 0.05, ** P < 0.01, *** P < 0.001, **** P < 0.0001.

Upon downloading and integrating expression data along with clinical data, no aberrant samples were identified through the sample dendrogram, and most clinical features were documented (Figure 1A). Subsequently, the WGCNA network was constructed using the gene expression profile after filtering out low-expressed genes. The WGCNA algorithm was applied to construct co-expression networks and modules for 472 samples from TCGA. We selected the top 25% mutated genes from the TCGA-SKCM cohort and performed clustering analysis using the ‘WGCNA’ R package. With a soft threshold β set to 11, the scale independence of the topological network exceeded 0.85, and the average connectivity approached 0 (Figures 1B, C). Therefore, setting the soft threshold to β=11 satisfied the scale-free topology criterion, with R^2 = 0.88 for the TOM (Figure 1D). The dynamic tree cutting algorithm based on TOM was employed to cluster all selected genes, resulting in the division of the tree into 12 modules (Figure 1E), each labeled with distinct colors.

Subsequently, we summarized the co-expression of feature genes and calculated the correlation between feature genes and clinical characteristics. Pathological diagnosis typically includes information on pathological stages M, N, and T at the time of diagnosis. Each module consists of distinct gene clusters, marked in an overlapping heatmap with different colors; red indicating positive correlation and green indicating negative correlation (Figure 2A). Among the co-expression modules and features we analyzed, we observed that the yellow module is most closely related to pathological stage M (R=0.43, P=1e-04); additionally, this module is positively correlated with age (R=0.25, P < 0.05) (Figure 2B). Therefore, we consider the yellow module containing 520 genes as a critical gene module and employ it for further analysis (Figure 2C).

To elucidate the functional mechanisms of the yellow module, we conducted enrichment analyses of GO terms and KEGG pathways. In this study, we initially examined biological processes (BP), molecular functions (MF), and cellular components (CC) to identify the most enriched GO pathways. The results revealed that, in terms of biological processes, these genes are primarily involved in epidermal development, keratinocyte differentiation, and keratinization processes. Regarding cellular components, the main enrichment was observed in cornified envelope, extracellular region, and extracellular region part. Molecular function analysis indicated that these genes are predominantly enriched in epidermal structural components (Figure 3A). Subsequently, based on the enrichment results of KEGG pathways, we identified pathways closely associated with tumorigenesis, including central carbon metabolism, cell adhesion molecules, junctional protein binding, and melanogenesis, providing important clues for understanding the functional mechanisms of genes within the yellow module (Figures 3B, C).

To explore the interrelations among genes and identify central genes, an initial analysis was conducted on genes within the yellow module, followed by visualization using STRING. After eliminating unconnected nodes, a PPI network was obtained (Figure 4A). Subsequently, employing ‘Cytoscape’, topological parameters for all nodes were computed, and the six highest-ranking node genes were identified as pivotal for subsequent analyses (Figure 4B). To investigate the association between central genes (DSC2, DSC3, DSG1, KRT6B, PKP1, and PKP3) and the response to oxidative stress, we employed GSVA to compute the oxidative stress scores for each sample. Pearson correlation analysis revealed a significant correlation between oxidative stress scores and DSC2 (R=0.45; P < 0.001) as well as DSC3 (R=0.094; P < 0.05), whereas other genes showed no involvement in the oxidative stress response (Figure 4C). This underscores the need for further exploration into the relationship between oxidative stress and DSC2, DSC3.

Through correlation analysis, we identified a set of 12 genes highly correlated with both DSC2 and DSC3, and delved into the interrelationships among them. Circular diagrams were presented for DSC2 (Figure 4D) and DSC3 (Figure 4E) separately to illustrate their associations. KM survival analysis was conducted on these selected central genes, revealing a significant association with the prognosis of SKCM patients (P < 0.05). As the expression levels of DSC2, DSC3, DSG1, KRT6B, PKP1, and PKP3 increased, the overall survival time of trunk subtype melanoma patients significantly decreased (Figure 5).

Immunohistochemical analysis revealed heightened infiltration levels of NK cells, with significant differences observed between groups with varying DSC2 expression levels. Likewise, significant disparities in dendritic cell activation and infiltration were evident between the two DSC3 expression level groups (Figure 6A). Furthermore, a positive correlation was identified between DSC2 expression and activation of CD4 memory T cells, dendritic cells, mast cells, and neutrophil infiltration. Conversely, DSC2 expression exhibited a negative correlation with regulatory T cells and activated NK cell infiltration (Figure 6B). Additionally, DSC3 expression showed a positive correlation with dendritic cell activation and neutrophil infiltration. Figure 6C depicts the association between DSC2 and DSC3 and highly correlated immune checkpoint genes.

Subsequently, we compared the functional pathways associated with different oxidative stress scores in melanoma. Utilizing GSVA, we assessed the oxidative stress-related signaling pathways in melanoma and further explored them based on KEGG and GO datasets. Enriched genes, signaling pathways, and functions were highlighted using heatmap plots. In the KEGG analysis, we investigated the top 20 associated signaling pathways such as the MAPK signaling pathway, JAK STAT signaling pathway (Figure 7A). Additionally, GO analysis revealed significant differential functional pathways between the two oxidative stress groups, presenting the top 20 signals in Figure 7B, including muscle cell proliferation, growth factor receptor binding.

We identified a close association between DSC2 and DSC3 with oxidative stress processes. Further investigation demonstrated a positive correlation between the expression of DSC2 and DSC3 and the infiltration of activated dendritic cells and neutrophils, while showing a negative correlation with memory B cells (Figures 8A, B). Subsequently, we delved deeper into the expression patterns of DSC2 and DSC3 at the single-cell level (Figure 8C). Remarkably, the average expression level of DSC2 was highest in monocytes, followed by pDCs (Figures 8D, E). Interestingly, DSC3 was predominantly enriched in NK cells, CD8 T cells, and CD4 T cells (Figures 8F, G), indicating distinct spatial distributions of these two oxidative stress molecules.

Immunotherapy, as a prominently featured anticancer strategy in recent years, has demonstrated significant efficacy in the treatment of SKCM. However, some patients exhibit insensitivity or even develop resistance to immunotherapy. Thus, we further investigated the role of oxidative stress-related genes in SKCM immunotherapy. It was observed that patients with high expression levels of DSC2 and DSC3 were more sensitive to immunotherapy (Figures 9A, B), possibly due to infiltration by macrophages and CD8 T cells (Figures 9C, D). Additionally, DSC2 could serve as an effective predictive biomarker for immunotherapy response in SKCM patients (Figure 10A). Intriguingly, SKCM patients with high expression of DSC2 exhibited poorer prognosis following Anti-PD-1/PD-L1 therapy compared to those with lower DSC2 expression (Figure 10B), which may be attributed to the pro-tumorigenic characteristics of DSC2. Furthermore, chemotherapy, as a crucial component of anticancer therapy, often fails ultimately due to the development of resistance. Therefore, we further investigated the expression of DSC2 in drug resistance. Results indicated that SKCM patients with high expression of DSC2 were more prone to developing resistance to Apitolisib, Motesanib, and Amuvatinib (Figure 11), suggesting DSC2 as one of the targeted strategies to enhance chemotherapy sensitivity.

Given the lack of expression profiles for adjacent normal tissues within the TCGA-SKCM cohort, a total of 92 SKCM samples were obtained, comprising both tumor and adjacent normal tissues. Differential expression analysis of these core proteins was subsequently performed through immunohistochemistry experiments (Figures 12A, 13A). The results indicate that in trunk subtype melanomas, the protein expression levels of DSC2, DSC3, DSG1, KRT6B, PKP1, and PKP3 are significantly higher than those in adjacent normal samples (Figures 12B, C), suggesting the crucial roles of these central genes in SKCM progression. However, in other subtype melanoma patients, only DSG1 and PKP1 exhibit significantly elevated protein expression levels compared to adjacent normal samples, with no observed differences in other genes between tumor and adjacent normal samples (Figures 13B, C).