Single-cell RNA analysis to identify five cytokines signaling in immune-related genes for melanoma survival prognosis

Melanoma is one of the deadliest skin cancers. Recently, developed single-cell sequencing has revealed fresh insights into melanoma. Cytokine signaling in the immune system is crucial for tumor development in melanoma. To evaluate melanoma patient diagnosis and treatment, the prediction value of cytokine signaling in immune-related genes (CSIRGs) is needed. In this study, the machine learning method of least absolute selection and shrinkage operator (LASSO) regression was used to establish a CSIRG prognostic signature of melanoma at the single-cell level. We discovered a 5-CSIRG signature that was substantially related to the overall survival of melanoma patients. We also constructed a nomogram that combined CSIRGs and clinical features. Overall survival of melanoma patients can be consistently predicted with good performance as well as accuracy by both the 5-CSIRG signature and nomograms. We compared the melanoma patients in the CSIRG high- and low-risk groups in terms of tumor mutation burden, infiltration of the immune system, and gene enrichment. High CSIRG-risk patients had a lower tumor mutational burden than low CSIRG-risk patients. The CSIRG high-risk patients had a higher infiltration of monocytes. Signaling pathways including oxidative phosphorylation, DNA replication, and aminoacyl tRNA biosynthesis were enriched in the high-risk group. For the first time, we constructed and validated a machine-learning model by single-cell RNA-sequencing datasets that have the potential to be a novel treatment target and might serve as a prognostic biomarker panel for melanoma. The 5-CSIRG signature may assist in predicting melanoma patient prognosis, biological characteristics, and appropriate therapy.

1 Introduction

Melanoma is the most dangerous type of skin cancer, accounting for 90% of all skin cancer deaths. Melanoma has become more common in recent decades, with an estimated 232,100 new cases and 55,500 deaths per year (1). Treatment of melanoma includes surgery, chemotherapy, radiotherapy, immunotherapy, targeted therapy, and other methods (2, 3). With a 90% cure rate, melanoma surgery remains the most effective treatment option. However, patients have a high rate of recurrence despite aggressive interventions, which contributes to the poor prognosis of melanoma (4). Therefore, new prognostic biomarkers for melanoma should be investigated to identify high-risk subpopulations and guide more effective individual treatments.

Lymphocytes, monocytes, macrophages, B cells, and T cells are immune system cells that create cytokines (5). Immune system cells can communicate with one another using cytokines to produce coordinated, efficient, but self-restraining antigen responses. Although there are numerous ways that the immune system can communicate with one another directly between cells, cytokine synthesis allows for a more varied and effective transmission of immunological information (6). Numerous malignancies, such as melanoma and renal cell carcinoma, are treated with cytokines (5). Cytokines at the tumor site stimulate immune effector cells, improving tumor cell recognition. The interaction between interleukins and interferons has an enhanced immunostimulatory impact (7). As a result, many cytokine- or cytokine antagonist-based cancer therapies have been developed (5, 8, 9). However, for patients with advanced-stage illness, cytokine-based therapy had a modest therapeutic effect.

Based on single-cell data and TCGA data, we used least absolute selection and shrinkage operator regression and Cox regression to build a prediction model of melanoma. The predictive significance of our cytokine signaling in immune-related genes (CSIRG) signature was further verified by receiver-operating characteristic analysis. Gene set enrichment analysis was conducted to help elucidate the intrinsic mechanisms. Moreover, the predictive significance of our CSIRG signature was confirmed using independent GEO data. Our results imply that the CSIRG signature plays a crucial role in predicting the prognosis of melanoma patients.

2 Methods

2.1 Data source

Melanoma single-cell datasets were accessed from GSE115978 (10). The TCGA (https://portal.gdc.cancer.gov/) (11) and GSE65904 cohorts (12), with 470 and 214 patients, respectively, provided bulk mRNA expression data in this study. Following data collection, an integrated analysis was conducted.

2.2 Data processing with single-cell RNA-sequencing

Seurat (version 4.1.0) was used to analyze the single-cell data of 31 melanoma patients. This cohort included 3 patients with primary melanoma and 28 patients with metastatic melanoma. Among these patients, 15 were untreated, and 16 were post-immunotherapy. The raw data were accessed from GSE115978 (10). Cells with > 600 and < 4,500 RNA features were included in the following analysis. Additionally, a linear dimensionality reduction was created using tSNE, and significant dimensions were identified with an estimated P value. Adjusted P < 0.05 and |logFC| > 0.25 were used for differentially expressed CSIRG analysis between metastatic melanoma and the primary melanoma group. We obtained cytokine signaling in immune-related genes (CSIRGs) from the Reactome pathway database.

2.3 Establishment of prognostic signature

The research included CSIRG signature genes with matching clinical characteristics. A training group (TCGA-SKCM) and a testing group (GSE65904) were formed from the melanoma patients. Univariate Cox regression was utilized in the TCGA-SKCM group to narrow down the prognostic CSIRG genes. The results from the univariate Cox regression were used for the least absolute selection and shrinkage operator (LASSO) regression. We further used the “glmnet” R package to perform the LASSO algorithm and select the potential candidates. A subsample of 1000 iterations was conducted on the dataset, and CSIRGs with occurrence frequencies of more than 950 were identified for further analysis. To narrow the list of candidate CSIRGs related to overall survival and construct a predictive gene signature, multivariate Cox regression was conducted. Finally, a risk score (RS) model of the CSIRG signature was constructed using the multivariate Cox results. Multivariate Cox regression was performed to acquire the regression coefficient (β), and the RS was generated using the coefficients and expressions. Using the findings, the RS equation is RS = coefficient1 * gene1 expression + coefficient2 * gene2 expression + coefficientN * geneN expression. Patients in the TCGA-SKCM and GSE65904 cohorts were divided into high- and low-risk groups by median RS. The survival analysis of melanoma patients was calculated using the Kaplan−Meier method. Receiver-operating characteristic (ROC) analyses validated the signature performance.

2.4 Survival analysis of 11 CSIRGs accessed from the results of LASSO regression

Kaplan–Meier survival analysis with the log-rank test was performed to study the 11 CSIRGs (ATF2, CCR1, CRKL, EIF4A2, IFI30, MCL1, NUP188, STAT1, STAT3, TNFSF13B, and YWHAZ) of melanoma patients in the TCGA-SKCM cohort. The patients were categorized into two high and low gene expression groups using the median expression level of each gene.

2.5 Nomogram construction

Before producing the nomogram, clinical and CSIRG signatures were combined. The appropriate clinical characteristics (including age, sex, tumor stage, tumor T stage, tumor M stage, and tumor N stage) and CSIRG RS were subsequently chosen using univariate and multivariate Cox regression models. The prognostic nomogram model was further developed using the CSIRG RS and independent clinical factors. ROC curves and decision curves were used to evaluate the performance of the nomogram model.

2.6 Mutation identification and tumor mutation burden quantification

TMB is base deletion, insertion, or substitution divided by the total number of variants divided by exon length. Maftools was used to analyze somatic mutation data to study TMB landscapes. The data were separated into two categories based on the CSIRG risk assessment. Mutant genes in the CSIRG high- and low-risk groups were analyzed.

2.7 Gene set enrichment analysis

Broad Institute’s software was used to conduct gene set enrichment analysis in the TCGA-SKCM cohort (13). A KEGG gene set (KEGG C2, MsigDB database v7.5.1) was used in this study. Gene set enrichment analysis was performed with 1,000 permutations (14). KEGG pathway study comparing the CSIRG high- and low-risk groups.

2.8 Immune infiltration analysis by xCell and CIBERSORT algorithm

We compared immune cell ratings using the xCell and CIBERSORT algorithms in the TCGA dataset (15, 16). To calculate the tumor-infiltrating cell distribution scores using xCell, this research integrated the transcriptome data of patients with the expression of marker genes from 64 different kinds of immune cells. Subgroup investigation of immune cell infiltration in the CSIRG high- and low-risk groups. To validate the xCell results, we further performed the most widely used immune infiltration algorithm of CIBERSORT to analyze the enrichment scores of 22 kinds of immune cells for each melanoma patient.

2.9 Statistical analysis

Statistical analyses and graphs were calculated using R version 4.0.5 and the necessary packages. P < 0.05 was considered to be statistically significant.

3 Results

3.1 Single-cell RNA sequencing analysis

We accessed the single-cell data of 31 patients from GSE115978 (10). Then, we performed single-cell transcriptome analysis to compare primary and metastatic melanoma by the Seurat pipeline. After quality control, a differential analysis was performed on the cells from the primary and metastatic melanoma tissues. The distribution of the single-cell data is shown in Figures 1A–H. The samples included 27 clusters (Figure 1A). We also showed the distribution of cells in different patient groups (Figure 1B, metastatic melanoma and primary melanoma groups). Malignant cells are categorized by patients (Figure 1C), while nonmalignant cells, including immune cells and stromal cells, are classified by their cell type (Figure 1D), as the original article reported. The nonmalignant cells included CD4+ T cells, CD8+ T cells, T cells, NK cells, macrophages cells, B cells, cancer-associated fibroblast (CAF) cells, and endothelial cells (Figure 1D). We further showed the distribution of malignant cells and nonmalignant cells in the posttreatment group (Figures 1E–F) and untreated group (Figures 1G–H). Among 31 patients, cells from only 23 patients were annotated with malignant cells which were shown in Figures 1C, E, F.

FIGURE 1

Figure 1 Single-cell RNA sequencing analysis of the GSE115978 dataset. (A) T-SNE of all of the melanoma samples with 27 clusters. (B) T-SNE of the metastatic melanoma and primary melanoma groups. (C) Malignant cells are categorized by the patients. (D) Nonmalignant cells, including immune cells and stromal cells, are classified by their cell type. The distribution of malignant cells and nonmalignant cells in the posttreatment group (E, F) and untreated group (G, H).

3.2 A risk model based on five cytokine signaling in immune-related genes

The CSIRGs used in the following analysis were obtained from the Reactome pathway database. As a consequence, 54 CSIRGs were differentially expressed between the metastatic melanoma and primary melanoma groups (Figure 2A). Compared with the primary melanoma group, 34 of the CSIRGs were upregulated in the metastatic melanoma group (Figure 2B; Supplementary Table 1), while 20 of the CSIRGs were downregulated in the metastatic melanoma group (Figure 2C; Supplementary Table 1). Among the upregulated CSIRGs, most were expressed in all of the cell types, including malignant cells (Mal), CD4+ T cells, CD8+ T cells, T cells, NK cells, macrophages cells, B cells, cancer-associated fibroblast (CAF) cells, and endothelial cells. However, TNFSF8 was highly expressed in CD4+ T cells and T cells of both the primary and metastatic groups. TNFSF4 was highly expressed in malignant cells and CD8+ T cells of the metastatic group. S100B and TRIM2 were highly expressed in malignant cells and CAF cells of the metastatic group. FN1 was highly expressed in Mal cells, CAF cells, and endothelial cells of the metastatic group. COL1A2 was highly expressed in CAF cells of the metastatic group. EGR1 was highly expressed in macrophage cells of the primary group. SHC1 was highly expressed in endothelial cells of the primary group (Figures 2B, C).

FIGURE 2

Figure 2 Identification of differentially expressed genes (DEGs) and cytokine signaling in immune-related genes (CSIRGs) between the metastatic melanoma and primary melanoma groups. (A) DEGs and CSIRGs in GSE115978 (|logFC| > 0.25 and adjusted P value < 0.05). (B) Dot plot showing the DEGs that were upregulated in the metastatic melanoma group. Mal: Malignant cells. (C) Dot plot showing the DEGs that were downregulated in the metastatic melanoma group. Mal: Malignant cells. (D, E) LASSO regression was performed on CSIRGs.

We performed univariate Cox regression analysis to screen prognostic CSIRGs from 54 CSIRGs. As a result, 20 CSIRGs were revealed to be highly linked with overall survival (OS) (Supplementary Table 2). LASSO regression further selected 11 CSIRGs (Figures 2D, E). Kaplan–Meier survival analysis of 11 CSIRGs, including ATF2, CCR1, CRKL, EIF4A2, IFI30, MCL1, NUP188, STAT1, STAT3, TNFSF13B, and YWHAZ, in melanoma patients in TCGA was performed. The results showed that high expression of CCR1, MCL1, STAT1, and TNFSF13B was correlated with longer survival time (Figure 3). High expression of NUP188 was correlated with shorter survival time (Figure 3). Then, multivariate Cox regression analysis was carried out to screen prognostic CSIRGs from 11 CSIRGs. Finally, five hub CSIRGs (EIF4A2, MCL1, NUP188, STAT1, and YWHAZ) were identified with the minimum Akaike information criterion (AIC) and were further used for CSIRG model construction (Table 1; Supplementary Table 3). The gene model established based on five hub genes was as follows: risk score (RS) = -0.00536× EIF4A2 - 0.00378 × MCL1 + 0.05332 × NUP188 - 0.00502 × STAT1 + 0.00571 × YWHAZ. The RS medians (TCGA-SKCM: 1.071; GSE65904: 1.023) were utilized to divide the TCGA-SKCM and GSE65904 datasets into two groups with high and low RS (Supplementary Tables 4, 5). The area under the receiver-operating characteristic curves (AUCs) proved that the gene model had relatively good accuracy. Specifically, the TCGA-SKCM dataset had AUCs of 0.693 (1 year), 0.669 (3 years), and 0.747 (5 years) (Figure 4A), whereas the GSE65904 dataset had AUCs of 0.669 (1 year), 0.645 (3 years), and 0.655 (5 years) (Figure 4B). Furthermore, Kaplan−Meier curves of both the TCGA-SKCM (Figure 4C) and GSE65904 datasets (Figure 4D) revealed that melanoma patients in the CSIRG high-risk group were associated with worse outcomes. The RS of each patient in the TCGA-SKCM cohort is displayed in Figure 4E. Figure 4F also confirmed that the mortality of the CSIRG high-risk melanoma patients increased. The expression of the five hub CSIRGs (EIF4A2, MCL1, NUP188, STAT1, and YWHAZ) in the high- and low-risk groups is shown in Figure 4G. The results showed that EIF4A2, MCL1, and STAT1 were downregulated in the high-risk group, while NUP188 and YWHAZ were upregulated in the high-risk group (Figure 4G).

FIGURE 3

Figure 3 Kaplan–Meier survival analysis of 11 CSIRGs, including ATF2, CCR1, CRKL, EIF4A2, IFI30, MCL1, NUP188, STAT1, STAT3, TNFSF13B, and YWHAZ, identified by LASSO regression in the TCGA-SKCM cohort.

TABLE 1

Table 1 Multivariate cox regression analysis of different combination of CSIRGs by the Akaike information criterion (AIC).

FIGURE 4

Figure 4 The prognostic value of the risk signature was validated. ROC curves for predicting 1-, 3-, and 5-year overall survival (OS) in the (A) TCGA and (B) GSE65904 datasets using the risk model. Kaplan−Meier curves were used to assess the OS probabilities in the training (C) and testing groups (D). (E) In the training dataset of TCGA, the risk score (RS) of each patient is displayed. (F) OS and survival status are shown in the TCGA cohort (red dots indicate death, and green dots indicate survival). (G) The expression of the five CSIRGs (EIF4A2, MCL1, NUP188, STAT1, and YWHAZ) used to construct the prediction model was analyzed in the high- and low-risk groups. : ***: P≤0.001 and ****: P≤0.0001.

3.3 Nomogram model construction for melanoma patients

To further construct a nomogram combining clinical variables (including age, sex, tumor stage, tumor T stage, tumor M stage, and tumor N stage) and CSIRG RS, univariate and multivariate Cox regression models were constructed from the TCGA-SKCM dataset (Table 2). As a result, the tumor T stage was selected and integrated with the RS to create a clinical nomogram using the generalized linear model regression technique (Figure 5A). The survival rates of melanoma patients over 1, 3, and 5 years could be predicted using the clinical nomogram. The AUCs for 1, 3, and 5 years were 0.81, 0.734, and 0.741 for the TCGA-SKCM dataset, respectively (Figure 5B). Decision curve analysis was carried out to evaluate the nomogram performance for 5-year OS (Figure 5C).

TABLE 2

Table 2 Univariate and multivariate Cox regression to analysis independent prognosis factors in the TCGA-SKCM cohort.

FIGURE 5

Figure 5 Development and evaluation of a risk-based predictive nomogram for melanoma. (A) A nomogram for predicting melanoma prognosis that includes the CSIRG RS and tumor T stage. (B) Receiver operating characteristic curves and the area under the receiver operating characteristic curve evaluation of the nomogram performance of 1-, 3-, and 5-year OS. (C) Decision curve evaluation of the nomogram performance for 5-year OS.

3.4 Tumor mutation burden analysis

The TMB profile of the TCGA-SKCM dataset was downloaded and matched with the RS for subsequent analysis. TTN, MUC16, DNAH5, BRAF, and PLCO were the top five mutated genes in the CSIRG high- and low-risk groups, and their mutation rates were higher in the low-risk group than in the high-risk group (Figures 6A, B). The sixth mutated gene in the CSIRG high-risk group was LRP1B, and the mutation rate of this gene was higher in the high-risk group (38%) than in the low-risk group (39%) (Figures 6A, B). Since TTN demonstrated the highest mutation (high risk: 69%, low risk: 74%) in the TCGA-SKCM dataset, it might be an important risk factor in melanoma (Figures 6A, B).

FIGURE 6

Figure 6 The difference in tumor mutation burden between the CSIRG low- and high-score groups is shown. (A, B) Gene mutation information in the high-risk (A) and low-risk (B) categories of the TCGA dataset.

3.5 Gene set enrichment analysis

To gain further insight into signaling pathways associated with melanoma, GSEA was performed on the TCGA-SKCM dataset (Supplementary Table 6). The results from GSEA revealed that aminoacyl tRNA biosynthesis, DNA replication, glyoxylate and dicarboxylate metabolism, oxidative phosphorylation, pyrimidine metabolism, and Vibrio cholerae infection were mainly enriched in the high-risk group (Figures 7A). In contrast, antigen processing and presentation, Leishmania infection, systemic lupus erythematosus, Toll-like receptor signaling pathway, type 1 diabetes mellitus, and viral myocarditis were mainly enriched in the low-risk group (Figures 7B).

FIGURE 7

Figure 7 KEGG and hallmark gene set enrichment analysis for the CSIRG high-risk (A) and low-risk (B) melanoma patient groups in the TCGA cohort.

3.6 Immune cell infiltration analysis

We then evaluated the influence of the RS on the immune microenvironment by comparing immune cell infiltration in the TCGA-SKCM dataset through the xCell algorithm. Increased enrichment of activated myeloid dendritic cells, B cells, memory CD4+ T cells, naïve CD4+ T cells, nonregulatory CD4+ T cells, naïve CD8+ T cells, CD8+ T cells, central memory CD8+ T cells, effector memory CD8+ T cells, class-switched memory B cells, common lymphoid progenitors, myeloid dendritic cells, macrophages, M1 macrophages, memory B cells, monocytes, naïve B cells, plasmacytoid dendritic cells, plasma B cells, gamma-delta T cells, Th2 CD4+ T cells, and regulatory T cells was found in the high-risk group, while increased enrichment of NKT cells and Th1 CD4+ T cells was found in the low-risk group (Figure 8A; Supplementary Table 7).

FIGURE 8

Figure 8 Correlations between melanoma RS and tumor immune microenvironment. The proportionate differences in immune and stromal cells between the CSIRG high- and low-risk melanoma patient groups were visualized using a violin plot by the xCell (A) and CIBERSORT (B) algorithms.

To validate the results of xCell, we further used the CIBERSORT algorithm to analyze immune cell infiltration. The results showed that increased enrichment of monocytes and M2 macrophages was found in the high-risk group, while increased enrichment of plasma cells was found in the low-risk group (Figure 8B). As a result, only monocytes were found to be positively correlated with the high-risk group in both the xCell and CIBERSORT algorithms.

4 Discussion

Previous research has shown that the prognoses of melanoma patients cannot be adequately predicted using routinely utilized clinicopathological markers. In this study, five cytokine signaling pathways in immune-related genes (CSIRGs, including EIF4A2, MCL1, NUP188, STAT1, and YWHAZ) were used to generate a predictive signature. This signature was used to classify the melanoma patients in the TCGA and GSE65904 cohorts into low-risk and high-risk groups, depending on their risk scores. In addition, the results of the receiver-operating characteristic (ROC) analysis demonstrated that the five-gene signature has good potential for predicting overall survival in melanoma patients. The expression of EIF4A2 has also been shown to have a favorable correlation with the prognosis of non-small cell lung cancer and breast cancer in a number of studies (17, 18). On the other hand, a recent study found that higher EIF4A2 expression in colorectal cancer (CRC) was linked to a worse chance of survival. In addition, the findings of research conducted on cells and animals have provided additional evidence that EIF4A2 has a role in both the promotion of CRC metastasis and oxaliplatin resistance (19). An antiapoptotic protein that promotes resistance to numerous chemotherapeutic agents (20) is encoded by the MCL1 gene, which is typically increased in melanoma as well as in many other kinds of tumors. Nucleoporin is a component of the nuclear pore complex, which is encoded by the NUP188 gene (NPC). In mitotic cells, nucleoporin plays a role in the regulation of chromosomal segregation by increasing chromosome alignment. It is possible that problems with the chromosomal segregation process are responsible for the aneuploidy that occurs in some cancer cells (21). Therefore, NUP188 might play a part in the process of oncogenesis. Studies have also demonstrated that the expression of STAT1 is lower in gliomas than in normal brain tissues (22, 23). This difference in expression has been proven to occur in gliomas. The overexpression of STAT1 significantly inhibits the development of glioma cells and stimulates apoptosis (24) in these cells. However, the function of aberrant STAT1 and the mechanism by which it causes melanoma are not yet fully understood. Emerging data have shown that YWHAZ is essential in the growth of many different kinds of tumors (25). YWHAZ has also been shown to be a useful prognostic marker for multiple cancers according to a number of studies (26–29).

We further analyzed the correlation of tumor mutation burden (TMB) in the CSIRG high- and low-risk groups. Patients in the high-TMB group had significantly better survival outcomes. The mutation of TTN, which can aid carcinogenesis and metastasis in melanoma, was proven to be the most common mutation in melanoma patients. We showed that TTN mutation might lead to a good overall outcome for melanoma patients in the low-risk group. The mutation of TTN, which should be investigated in further studies, was negatively correlated with RS in our study.

The association of immune cell infiltration with risk score was explored, and the modulation of immune cells by CSIRGs is relatively underexplored in skin melanoma. We compared the immunological microenvironments (TMEs) of the CSIRG high- and low-risk groups and discovered that the former was enriched with monocytes. Although monocytes were more enriched in the high-risk group, survival was worse in the high-risk group. Previous studies proved that there are functionally distinct subsets of monocytes in different tumor microenvironments (30, 31). Pathologically activated monocytes are myeloid-derived suppressor cells (MDSCs) that have potent immunosuppressive effects (32). Important differences between MDSCs and classical monocytes were also reported (32, 33). In addition to causing tumor vasculogenesis, MDSCs suppress T-cell immunity and enable tumors to escape immunity (33). (PMID: 24060865). A previous study also proved that a poor prognosis was associated with phosphorylated STAT3 expression in monocytes (34). Aberrantly hyperactivated STAT3 monocytes promote liver tumorigenesis in both clinical patients and in vivo animal experiments (34). It is possible that the enriched monocytes in the CSIRG high-risk group have the same characteristics as MDSCs, which have potent immunosuppressive activity.

CSIRG may participate in the TME and immunological responses. Among the five CSIRGs, STAT1 and MCL1 were previously reported to be associated with infiltrating immune cells. JAK1/STAT1 activation occurs in immune-exhausted cells and is associated with increased Treg cell scores and upregulation of PDL1 (35). STAT1 in the majority of immune cell populations was more active in the high-risk group, but survival was worse. It was also reported that STAT1-deficient mice exhibited decreased accumulation of Th1 cells (36). Moreover, STAT1 is a biomarker of immune infiltration changes after anti-tuberculosis treatment (37). STAT1 was positively correlated with monocytes and neutrophils and negatively correlated with CD8+ T cells (37). In a study of acute lung injury, MCL1 was found to be negatively correlated with both B cells and T cells (38). Further in-depth research is required to study the relationship between CSIRG and the TME.

We further applied GSEA to detect the genet set enrichment in the CSIRG high- and low-risk groups. KEGG pathway analysis showed that these genes are associated with cytokine signaling in immune pathways. Enrichment with phenotypic consistency was also found in pathways such as glyoxylate and dicarboxylate metabolism, oxidative phosphorylation, DNA replication, aminoacyl tRNA biosynthesis, and pyrimidine metabolism, which were enriched in the high-risk group, whereas antigen processing and presentation, the Toll-like receptor signaling pathway and viral myocarditis were enriched in the low-risk group. Above all, it is reasonable to suggest that these CSIRGs might participate in the occurrence and development of melanoma through these pathways.

In summary, we defined a novel CSIRG signature in melanoma. This CSIRG signature and nomogram showed high predictive capability for the prognosis of melanoma patients. New insights into the prognosis of patients with melanoma may be provided by the CSIRG signature. Melanoma therapeutic targets may be developed through pathways related to CSIRG. Patients may benefit from further screening of anticancer drugs sensitive to both high and low CSIRG groups.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Author contributions

ZP and LM initiated the study. QZ and LM performed the data analysis and manuscript writing. JC and YX revised the manuscript. XZ supervised the study. All authors contributed to the article and approved the submitted version.

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

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

Abbreviations

AIC, Akaike information criterion; AUCs, The area under the receiver operating characteristic curves; CSIRGs, Cytokine signaling in immune-related genes; DEGs, Differentially expressed genes; LASSO, Least absolute selection and shrinkage operator; OS, Overall survival; ROC, Receiver-operating characteristic; RS, Risk score; TMB, Somatic tumor mutational burden; GSEA, Gene set enrichment analysis.

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Davis LE, Shalin SC, Tackett AJ. Current state of melanoma diagnosis and treatment. Cancer Biol Ther (2019) 20:1366–79. doi: 10.1080/15384047.2019.1640032

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Michielin O, van Akkooi ACJ, Ascierto PA, Dummer R, Keilholz U, ESMO Guidelines Committee. Cutaneous melanoma: ESMO clinical practice guidelines for diagnosis, treatment and follow-up†. Ann Oncol (2019) 30:1884–901. doi: 10.1093/annonc/mdz411

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Eggermont AMM, Hamid O, Long GV, Luke JJ. Optimal systemic therapy for high-risk resectable melanoma. Nat Rev Clin Oncol (2022) 19:431–9. doi: 10.1038/s41571-022-00630-4

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Propper DJ, Balkwill FR. Harnessing cytokines and chemokines for cancer therapy. Nat Rev Clin Oncol (2022) 19:237–53. doi: 10.1038/s41571-021-00588-9

PubMed Abstract | CrossRef Full Text | Google Scholar

6. 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:541–50. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Lee S, Margolin K. Cytokines in cancer immunotherapy. Cancers (Basel) (2011) 3:3856–93. doi: 10.3390/cancers3043856

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Myers JA, Miller JS. Exploring the NK cell platform for cancer immunotherapy. Nat Rev Clin Oncol (2021) 18:85–100. doi: 10.1038/s41571-020-0426-7

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Saxton RA, Glassman CR, Garcia KC. Emerging principles of cytokine pharmacology and therapeutics. Nat Rev Drug Discovery (2023) 22:21–37. doi: 10.1038/s41573-022-00557-6

CrossRef Full Text | Google Scholar

10. Jerby-Arnon L, Shah P, Cuoco MS, Rodman C, Su M-J, Melms JC, et al. A cancer cell program promotes T cell exclusion and resistance to checkpoint blockade. Cell (2018) 175:984–997.e24. doi: 10.1016/j.cell.2018.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wheeler DA, Roberts LR. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell (2017) 169:1327–1341.e23. doi: 10.1016/j.cell.2017.05.046

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Cirenajwis H, Ekedahl H, Lauss M, Harbst K, Carneiro A, Enoksson J, et al. Molecular stratification of metastatic melanoma using gene expression profiling: Prediction of survival outcome and benefit from molecular targeted therapy. Oncotarget (2015) 6:12297–309. doi: 10.18632/oncotarget.3655

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (2011) 27:1739–40. doi: 10.1093/bioinformatics/btr260

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U.S.A. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Aran D, Hu Z, Butte AJ. XCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

16. 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:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Shaoyan X, Juanjuan Y, Yalan T, Ping H, Jianzhong L, Qinian W. Downregulation of EIF4A2 in non-small-cell lung cancer associates with poor prognosis. Clin Lung Cancer (2013) 14:658–65. doi: 10.1016/j.cllc.2013.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yan LX, Wu QN, Zhang Y, Li YY, Liao DZ, Hou JH, et al. Knockdown of miR-21 in human breast cancer cell lines inhibits proliferation, in vitro migration and in vivo tumor growth. Breast Cancer Res (2011) 13:R2. doi: 10.1186/bcr2803

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Chen Z-H, Qi J-J, Wu Q-N, Lu J-H, Liu Z-X, Wang Y, et al. Eukaryotic initiation factor 4A2 promotes experimental metastasis and oxaliplatin resistance in colorectal cancer. J Exp Clin Cancer Res (2019) 38:196. doi: 10.1186/s13046-019-1178-z

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Belmar J, Fesik SW. Small molecule mcl-1 inhibitors for the treatment of cancer. Pharmacol Ther (2015) 145:76–84. doi: 10.1016/j.pharmthera.2014.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Itoh G, Sugino S, Ikeda M, Mizuguchi M, Kanno S, Amin MA, et al. Nucleoporin Nup188 is required for chromosome alignment in mitosis. Cancer Sci (2013) 104:871–9. doi: 10.1111/cas.12159

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Chen J-J, Ren Y-L, Shu C-J, Zhang Y, Chen M-J, Xu J, et al. JP3, an antiangiogenic peptide, inhibits growth and metastasis of gastric cancer through TRIM25/SP1/MMP2 axis. J Exp Clin Cancer Res (2020) 39:118. doi: 10.1186/s13046-020-01617-8

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhang Y, Jin G, Zhang J, Mi R, Zhou Y, Fan W, et al. Overexpression of STAT1 suppresses angiogenesis under hypoxia by regulating VEGF−A in human glioma cells. BioMed Pharmacother (2018) 104:566–75. doi: 10.1016/j.biopha.2018.05.079

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Ehrmann J, Strakova N, Vrzalikova K, Hezova R, Kolar Z. Expression of STATs and their inhibitors SOCS and PIAS in brain tumors. In vitro and in vivo study. Neoplasma (2008) 55:482–7.

PubMed Abstract | Google Scholar

25. Chen C-H, Chuang S-M, Yang M-F, Liao J-W, Yu S-L, Chen JJW. A novel function of YWHAZ/β-catenin axis in promoting epithelial-mesenchymal transition and lung cancer metastasis. Mol Cancer Res (2012) 10:1319–31. doi: 10.1158/1541-7786.MCR-12-0189

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Lu J, Guo H, Treekitkarnmongkol W, Li P, Zhang J, Shi B, et al. 14-3-3zeta cooperates with ErbB2 to promote ductal carcinoma in situ progression to invasive breast cancer by inducing epithelial-mesenchymal transition. Cancer Cell (2009) 16:195–207. doi: 10.1016/j.ccr.2009.08.010

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Fan T, Li R, Todd NW, Qiu Q, Fang H-B, Wang H, et al. Up-regulation of 14-3-3zeta in lung cancer and its implication as prognostic and therapeutic target. Cancer Res (2007) 67:7901–6. doi: 10.1158/0008-5472.CAN-07-0090

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Matta A, DeSouza LV, Shukla NK, Gupta SD, Ralhan R, Siu KWM. Prognostic significance of head-and-neck cancer biomarkers previously discovered and identified using iTRAQ-labeling and multidimensional liquid chromatography-tandem mass spectrometry. J Proteome Res (2008) 7:2078–87. doi: 10.1021/pr7007797

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Nishimura Y, Komatsu S, Ichikawa D, Nagata H, Hirajima S, Takeshita H, et al. Overexpression of YWHAZ relates to tumor cell proliferation and malignant outcome of gastric carcinoma. Br J Cancer (2013) 108:1324–31. doi: 10.1038/bjc.2013.65

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kwiecień I, Rutkowska E, Raniszewska A, Rzepecki P, Domagała-Kulawik J. Modulation of the immune response by heterogeneous monocytes and dendritic cells in lung cancer. World J Clin Oncol (2021) 12:966–82. doi: 10.5306/wjco.v12.i11.966

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Ugel S, Canè S, De Sanctis F, Bronte V. Monocytes in the tumor microenvironment. Annu Rev Pathol (2021) 16:93–122. doi: 10.1146/annurev-pathmechdis-012418-013058

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Veglia F, Sanseviero E, Gabrilovich DI. Myeloid-derived suppressor cells in the era of increasing myeloid cell diversity. Nat Rev Immunol (2021) 21:485–98. doi: 10.1038/s41577-020-00490-y

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Talmadge JE, Gabrilovich DI. History of myeloid-derived suppressor cells. Nat Rev Cancer (2013) 13:739–52. doi: 10.1038/nrc3581

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wu W-Y, Li J, Wu Z-S, Zhang C-L, Meng X-L. STAT3 activation in monocytes accelerates liver cancer progression. BMC Cancer (2011) 11:506. doi: 10.1186/1471-2407-11-506

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Zhang F, Zhang Q, Zhu J, Yao B, Ma C, Qiao N, et al. Integrated proteogenomic characterization across major histological types of pituitary neuroendocrine tumors. Cell Res (2022) 32:1047–67. doi: 10.1038/s41422-022-00736-5

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Mikhak Z, Fleming CM, Medoff BD, Thomas SY, Tager AM, Campanella GS, et al. STAT1 in peripheral tissue differentially regulates homing of antigen-specific Th1 and Th2 cells. J Immunol (2006) 176:4959–67. doi: 10.4049/jimmunol.176.8.4959

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Zhou C, Liang T, Jiang J, Chen J, Chen T, Huang S, et al. MMP9 and STAT1 are biomarkers of the change in immune infiltration after anti-tuberculosis therapy, and the immune status can identify patients with spinal tuberculosis. Int Immunopharmacol (2023) 116:109588. doi: 10.1016/j.intimp.2022.109588

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ling D, Zhang X, Wu J, Xu Q, He Z, Zhang J. Identification of immune infiltration and effective immune biomarkers in acute lung injury by bioinformatics analysis. Cell Transplant (2022) 31:9636897221124484. doi: 10.1177/09636897221124485

CrossRef Full Text | Google Scholar