Sage Journals: Discover world-class research

Abstract

BACKGROUND:

Melanoma is fatal cancer originating from melanocytes, whose high metastatic potential leads to an extremely poor prognosis.

OBJECTIVE:

This study aimed to reveal the relationship among EMT, TIICs, and immune checkpoints in melanoma.

METHODS:

Gene expression data and clinical data of melanoma were downloaded from TCGA, UCSC Xena and GEO databases. EMT-related DEGs were detected for risk score calculation. “ESTIMATE” and “xCell” were used for estimating TIICs and obtaining 64 immune cell subtypes, respectively. Moreover, we evaluated the relationship between the risk score and immune cell subtypes and immune checkpoints.

RESULTS:

Seven EMT-related genes were selected to establish a risk scoring system because of their integrated prognostic relevance. The results of GSEA revealed that most of the gene sets focused on immune-related pathways in the low-risk score group. The risk score was significantly correlated with the xCell score of some TIICs, which significantly affected the prognosis of melanoma. Patients with a low-risk score may be associated with a better response to ICI therapy.

CONCLUSION:

The individualized risk score could effectively conduct risk stratification, overall survival prediction, ICI therapy prediction, and TME judgment for patients with melanoma, which would be conducive to patients’ precise treatment.

Keywords

EMT TIICs TME biomarkers immune checkpoints

Abbreviations

TCGA	The Cancer Genome Atlas
SKCM	skin cutaneous melanoma
EMT	endothelial-mesenchymal transition
GSEA	Gene Set Enrichment Analysis
GO	gene ontology
KEGG	Kyoto Encyclopedia of Genes and
	Genomes

DEGs	differentially expressed genes
STAT1	Signal Transducer and Activator of
	Transcription 1
ESRP1	Epithelial Splicing Regulatory
	Protein 1
PIM2	Pim-2 Proto-Oncogene
JAG2	Jagged Canonical Notch Ligand 2
FBP1	Fructose-Bisphosphatase 1
NME1	NM23 Nucleoside Diphosphate
	Kinase 1
PSME3	Proteasome Activator Subunit 3
TIICs	tumor infiltration immune cells
TME	Tumor microenvironment

ICI	immune checkpoint inhibitor
PD-1	programmed death 1
PD-L1	programmed death 1 ligand

1. Introduction

Melanoma is fatal cancer originating from melanocy- tes, which comprises 65% of all skin cancers and whose high metastatic potential leads to an extremely poor prognosis [1, 2]. Several factors such as ultraviolet, radiation exposure, and malignant moles transformation contribute to melanoma formation [3]. In addition, the expression of oncogenes is associated with melanoma [4].

A previous study reported that the reactivation of epithelial-mesenchymal transition (EMT) was associated with the infiltration and metastasis of melanoma [5]. EMT promotes melanoma metastatic potential through epithelial trait loss and mesenchymal characteristics acquisition. In addition, the tumor microenvironment (TME) that is associated with tumor growth and metastasis has a relationship with tumor prognosis, including melanoma. Moreover, a strong association had been reported between EMT status and tumor immune microenvironment in melanoma [6, 7]. EMT-related differentially expressed genes (DEGs) in the TME are involved in regulating the function of tumor infiltration immune cells (TIICs), which has direct killing effects on tumor cells that may be the reason of tumor progression and metastasis.

A recent study demonstrated that programmed cell death protein 1 (PD-1) therapy-resistant patients with melanoma present distinct signatures of upregulated genes involved in monocyte chemotaxis and EMT [8]. The extracranial activation/release of CD8 $+$ T cells from PD-1/CTLA-4 inhibition and potentiation of their recruitment to the brain are paramount to the intracranial anti-PD-1/CTLA-4 activity, suggesting enhancement of these processes as an immune therapy-enhancing strategy in metastatic brain cancer [9]. These data collectively highlight the role of EMT-related gene expression as a driver of phenotype switching in melanoma cells, which may affect the development of antitumor immunity and response to immune checkpoint blockade. Some studies have suggested that patients with tumors with higher EMT-related gene expression are more likely to benefit from immune checkpoint inhibitor (ICI) therapy [10]. However, how EMT affects the TME of melanoma is unclear.

In this study, we constructed a risk model on the basis of the Cancer Genome Atlas (TCGA) cohort and validated it using the GSE65904 cohort. The risk model is based on seven EMT-related genes (ESRP1, FBP1, JAG2, NME1, PIM2, PSME3, and STAT1) which were found to be strongly associated with CD8 $+$ T cells subtypes and dendritic cells (DCs) subtypes in the TME. Furthermore, the risk score was negatively related to PD-1, programmed cell death-ligand 1 (PD-L1), and CTLA-4, which were immune checkpoint proteins that affect the survival rate and efficacy of the immune checkpoints blockade. Our risk model can predict the efficacy of immune checkpoints blockade in melanoma. Furthermore, the nomogram based on the risk score and clinical characteristics could precisely predict the survival probability.

2. Patients and methods

2.1 Microarray data acquisition

We downloaded gene expression data and clinical data of skin cutaneous melanoma from the TCGA database (https://tcga-data.nci.nih.gov/tcga/). For further investigation of the underlying mechanisms in melanoma, fragments per kilobase of exon per million reads mapped reads data of melanoma in TCGA, and normal tissues in genotype tissue expression database were downloaded from the UCSC Xena database (https://xenabrowser.net/), which is free of computational batch effects. GSE65904 was downloaded from the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/gds). The EMT-related gene list was downloaded from the EMT gene database (http://dbemt.bioinfo-minzhao.org/download.cgi), which was used to collect information on human EMT-related genes [11].

2.2 Identification of DEGs and PPI network construction

“Estimate” package (version 1.0.13) R language (version 4.0.0) was used to estimate the ratio of the immune-stromal component in TME for each sample. TME results included the immune and stromal scores. The higher the respective score is, the larger the ratio of the corresponding component in TME would be. The EMT-related DEGs were obtained by using the R software “limma” package according to the parameters $|$ Log (fold change) $|$ $<$ 1 and adjusted $p$ -value $<$ 0.05. The common DEGs shared by melanoma, immune score, and metastatic status were selected for further analysis. DEGs were performed for enrichment analyses using the “clusterProfiler” package (version 3.16.0) and then visualized using the “ggplot2” package (version 3.3.0) in R studio. Gene oncology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) terms that adjusted $p$ -value $<$ 0.05 was considered significantly enriched. The protein-protein interaction (PPI) network was predicted using the search tool for the retrieval of interacting genes (http://string-db.org) (version 11.0) online database and visualized using Cytoscape (version 3.6.1).

2.3 Establishment of a prognostic gene signature for melanoma

Univariate Cox proportional hazards regression analysis was used to identify prognostic relevant genes from EMT-related DEGs. Moreover, Kaplan-Meier plots and log-rank tests were adopted to validate the relationship between the DEGs and overall survival. The 20 out of the 28 DEGs obtained from the Cox regression and Kaplan-Meier analyses were further used for LASSO regression analysis to select the genes included in the model. Risk scores were calculated using the generated coefficients and corresponding expression. Patients were divided into “high-risk” and “low-risk” groups on the basis of the median risk score obtained using this prognostic model. To further assess the diagnostic value of the risk score system, time-dependent receiver operating characteristics (ROC) and survival analysis were performed.

2.4 Gene set enrichment analysis (GSEA)

Hallmarks (h.all.v72.symbols.gmt) and c2 (c2.cp.kegg.v7.2.symbols.gmt) gene sets were downloaded from the Molecular Signatures Database [12] as target sets, using which gene set enrichment analysis (GSEA; V4.0.3) was performed between high-risk and low-risk score groups. The significant cutoff value was defined as the false discovery rate of $<$ 0.25 and the normal $p$ value $<$ 0.05.

2.5 xCell score estimation

The “xCell” package was used to calculate the cell type enrichment score (xCell scores) based on the gene expression data. Notably, the xCell tool provides 64 cell types, including lymphoid, myeloid, stromal cells, stem cells, and other cells. The difference between each cell subtype between the high- and low-risk score groups was compared using the $t$ -test method. The cutoff $p$ -value for the cell analyses was $<$ 0.05. The common cell types from the two datasets were compared using Venn diagrams.

2.6 Relationships between risk score, immune cells and immune checkpoints

Pearson’s correlation analysis between the risk score and the xCell score of immunity cells was performed using the “corrplot” package (version 0.84). Pearson’s correlation coefficient of $>$ 0.5 and $p$ -value of $<$ 0.05 were considered as significantly correlated. Moreover, Kaplan-Meier plots and log-rank tests were used to validate the relationship between the common immunity cells and overall survival. Furthermore, the difference between each immune checkpoint expression between the high- and low-risk score groups was compared using the $t$ -test method. The patients were divided into four groups on the basis of the risk score and immune checkpoints’ expression, and their relationships with overall survival were analyzed using Kaplan-Meier analysis.

2.7 Construction of nomograms

The nomograms were constructed using patients with melanoma in TCGA dataset. Moreover, univariate and multivariate Cox proportional hazards regression analyses were performed on the risk score of EMT-related genes signature and other clinical characteristics. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated. Nomograms were constructed in this study using information obtained from the results of multivariate Cox regression analysis. The predictive accuracy of the nomogram was assessed by ROC curve analysis.

2.8 Immunohistochemistry

Patients volunteering for the study had signed informed consent forms. In addition, this study was approved by the Ethics Department of the First Affiliated Hospital of Guangxi Medical University and conformed to the World Medical Association Declaration of Helsinki. Six pairs (melanoma and paracancerous tissues) of pathological sections for each gene used for immunohistochemistry analysis. The PSME3 and NME1 antibodies for immunohistochemical staining were purchased from the ABclonal. The PIM2, JAG2 and ESRP1antibodies were purchased from the Abmart. The FBP1 and STAT1 antibodies were purchased from the Abcam. After paraffin removal, hydration, and sealing, the specimens were mixed with anti-ESRP1, -FBP1, -JAG2, -NME1, -PIM2, -PSME3, and -STAT1 and incubated overnight at 4 ${}^{\circ}$ C (dilution ratios of 1:250, 1:200, 1:4000, 1:300, 1:100) [13]. We performed the statistical analysis of the positivity rate of the immunohistology images using Image J software and GraphPad Prism 8.

Table 1
Univariate Cox regression analyses to the DEGs

Genes	Full name	HR	CI95%	$P$ value
STAT1	Signal Transducer and Activator of Transcription 1	0.99	0.99–1	$P<$ 0.0001
HAVCR2	Hepatitis A Virus Cellular Receptor 2	0.94	0.91–0.97	6.00E-05
CCR5	C-C Motif Chemokine Receptor 5	0.93	0.9–0.96	8.00E-05
CXCL9	C-X-C Motif Chemokine Ligand 9	1	0.99–1	9.00E-05
TNFAIP8L2	TNF Alpha Induced Protein 8 Like 2	0.94	0.91–0.97	1.60E-04
CCL5	C-C Motif Chemokine Ligand 5	1	0.99–1	1.90E-04
PDCD1	Programmed Cell Death 1	0.94	0.91–0.97	2.00E-04
CD14	CD14 Molecule	0.99	0.99–1	2.90E-04
PIM2	Pim-2 Proto-Oncogene	0.99	0.98–0.99	4.50E-04
JAK3	Janus Kinase 3	0.95	0.92–0.98	7.30E-04
CXCL13	C-X-C Motif Chemokine Ligand 13	0.99	0.98–0.99	7.50E-04
LYN	LYN Proto-Oncogene	0.95	0.92–0.98	9.00E-04
CCR2	C-C Motif Chemokine Receptor 2	0.88	0.81–0.95	1.61E-03
VSIG4	V-Set and Immunoglobulin Domain Containing 4	0.98	0.97–0.99	2.07E-03
FBP1	Fructose-Bisphosphatase 1	0.97	0.95–0.99	2.08E-03
ESRP1	Epithelial Splicing Regulatory Protein 1	1.02	1.01–1.04	3.63E-03
IRF8	Interferon Regulatory Factor 8	0.96	0.94–0.99	5.01E-03
JAG2	Jagged Canonical Notch Ligand 2	1.05	1.01–1.09	6.16E-03
NME1	NME/NM23 Nucleoside Diphosphate Kinase 1	1.01	1–1.02	8.30E-03
PSME3	Proteasome Activator Subunit 3	1.02	1–1.03	1.31E-02
MMP9	Matrix Metallopeptidase 9	1	1–1	2.71E-02
UCP2	Uncoupling Protein 2	0.99	0.99–1	5.78E-02
KLF5	Kruppel Like Factor 5	1.03	1–1.05	6.00E-02
CCR7	C-C Motif Chemokine Receptor 7	0.98	0.97–1	7.66E-02
CXCR4	C-X-C Motif Chemokine Receptor 4	1	0.99–1	8.75E-02
CCL19	C-C Motif Chemokine Ligand 19	1	1–1	2.06E-01
CCL21	C-C Motif Chemokine Ligand 21	1	1–1	3.33E-01
CIRBP	Cold Inducible RNA Binding Protein	1	0.99–1	3.90E-01

3. Results

3.1 EMT-related DEGs were mainly involved in inflammation response

“The raw score obtained from “Estimate” package for the TCGA cohort was shown in supplementary file S1. Patients had a better prognosis in the high immune score group than the low immune score group (Fig. 1A). However, patients with a high stromal score had no effect on prolonging the prognosis (Fig. 1B). A total of 801 DEGs (168 upregulated and 633 downregulated) were identified in melanoma compared with the normal skin (Fig. 1C). In addition, 256 DEGs (135 upregulated and 121 downregulated) were identified in the high immune score group compared with the low immune score group (Fig. 1D). Moreover, 297 DEGs (228 upregulated and 69 downregulated) were identified in the metastatic melanoma compared with nonmetastatic melanoma (Fig. 1E). The common DEGs (22 up-regulated and 6 down-regulated) shared by melanoma, immune score, and metastatic status were selected for further analysis (Fig. 1F and G). The PPI network is shown in Fig. 1H. GO and KEGG enrichment analyses showed that DEGs were mainly involved in inflammation response, including T-cell activation and response to chemokine and cytokine (Supplementary Fig. S1A, S1B; Supplementary File S2).

Figure 1.

Correlation of Immune status and the prognosis in melanoma and DEGs identification. (A) Kaplan-Meier survival curve for the immune score with $p<$ 0.0001 by log-rank test. (B) Kaplan-Meier survival curve for the stromal score with $p=$ 0.02 by log-rank test. (C) Volcano plots for DEGs were generated by comparison of the melanoma samples vs. the normal skin samples. (D) Volcano plots for DEGs were generated by comparison of the high score group vs. the low score group in the immune score. (E) Volcano plots for DEGs were generated by comparison of the metastatic melanoma samples vs. primary melanoma samples. (F) Venn plots showing common up-regulated DEGs shared by melanoma, immune score and metastatic status. (G) Venn plots showing common down-regulated DEGs shared by melanoma, immune score and metastatic status. (H) Protein-protein interaction network of common DEGs in the TCGA. Red nodes represent up-regulated genes and green nodes represent down-regulated genes.

Figure 2.

Establishment of a prognostic gene signature for melanoma. (A, B) 7 of the DEGs selected by LASSO Cox regression analysis. Left: using 1000-fold cross-validation to the optimal penalty parameter lambda. Right: LASSO coefficient profiles of the 20 DEGs. (C) Distribution of the risk score in the TCGA cohort. Upper panel: classification of patients into different risk groups based on the median risk score. Middle panel: distribution of patients’ survival time and status. Bottom panel: Heatmap of expression profiles of included risk score related genes. (D) Kaplan-Meier survival curves between low- and high-risk groups. (E, F) Validation of the risk score system prognostic signature in the GSE65904.

Figure 3.

GSEA analysis. (A) GSEA in the Hallmarks (left) and KEGG (right) gene set between group high ( $n=$ 229) and low ( $n=$ 229) of the risk score in the TCGA cohort. (B) GSEA in the Hallmarks (left) and KEGG (right) gene set between group high ( $n=$ 105) and low ( $n=$ 105) of the risk score in the GSE65904 cohort.

3.2 A seven EMT-related genes signature was established

A total of 21 out of 28 DEGs were screened as the significant prognosis factors from the univariate Cox regression analysis (Table 1). Moreover, 24 out of 28 DEGs affected melanoma prognosis as noted from Kaplan-Meier analysis (Supplementary Fig. S2A). The 20 out of 28 DEGs shared by Cox regression analysis and Kaplan-Meier analysis (Supplementary Fig. S2B) were further used to perform LASSO regression analysis for selecting the genes included in the model (Fig. 2A and B). A seven EMT-related DEGs (ESRP1, FBP1, JAG2, NME1, PIM2, PSME3, and STAT1) risk signature was established. All these seven DEGs were significantly differently expressed between high- and low-risk score groups in the TCGA (Supplementary Fig. S3A) and GSE65904 cohort (Supplementary Fig. S3B). On the basis of the median risk score, patients were separated into two groups, high- and low-risk groups, wherein the high-risk group had lower survival rates (Fig. 2C and D). The GSE65904 cohort, which was used to validate the prognostic model, showed similar results (Fig. 2E and F).

3.3 The risk score predicts the involvement of immune pathways

A large number of genes were enriched in immunity- and tumor-related pathways. Results from the GSEA analysis were similar to GO and KEGG enrichment analyses to the DEGs. In the TCGA cohort, the Hallmark background showed enrichment mainly for apoptosis, IL2/STAT5 signaling, IL6/JAK/STAT3 signaling, inflammatory response, interferon $\alpha$ / $\gamma$ response, PI3K/AKT MTOR signaling, TGF beta signaling, and TNFa signaling through NF $\kappa$ B, whereas the KEGG background showed enrichment mainly for antigen processing and presentation, apoptosis, B-cell receptor signaling pathway, cell adhesion molecules, chemokine signaling pathway, cytokine-cytokine receptor interaction, JAK/STAT signaling pathway, MAPK signaling pathway, natural killer cell-mediated cytotoxicity, T-cell receptor signaling pathway, and Toll-like receptor signaling pathway (Fig. 3A). The mainly enriched pathways for Hallmark and KEGG backgrounds in the GSE65904 cohort (Fig. 3B) were in accordance with that in the TCGA cohort.

3.4 TIICs profile in different risk groups

The lymphoid and myeloid cells were considerably lower in the high-risk score group in the TCGA cohort (Supplementary Fig. S4A, S4B, S4C, S4D, S4E) and GSE65904 cohort (Supplementary Fig. S5A, S5B, S5C, S5D, S5E). Three subtypes of immunity cells were significantly upregulated, whereas 25 subtypes of immunity cells were significantly downregulated in the high-risk score group (Supplementary Fig. S5F). Pearson’s correlation analysis showed that 14 and eight subtypes of TIICs were significantly correlated with the risk score in the TCGA (Fig. 4A) and GSE65904 cohorts (Fig. 4B), respectively. The Venn plot showed that activated dendritic cells (aDC), pDC, CD8 $+$ T cells, CD8 $+$ Tcm, and CD8 $+$ Tem were considered TIICs in melanoma (Fig. 5A). Furthermore, the infiltration of aDC, pDC, CD8 $+$ T cells, and CD8 $+$ Tem showed a better prognosis in the TCGA (Fig. 5B–F) and GSE65904 cohorts (Fig. 6A–E) in the melanoma.

Figure 4.

The correlation of risk score with xCell score. (A) The correlation between the risk score and xCell score in the TCGA cohort. (B) The correlation between the risk score and xCell score in the GSE65904 cohort.

Figure 5.

Kaplan-Meier survival analysis to the immunity cells. (A) Venn plot showing the immunity cells that was significantly correlated with the risk score in the TCGA cohort and the GSE65904 cohort. Kaplan-Meier survival curves between low- and high-aDC(B), pDC(C), CD8 $+$ T-cells(D), CD8 $+$ Tcm(E), CD8 $+$ Tem(F) xCell score in the TCGA cohort.

Figure 6.

Kaplan-Meier survival curves between low- and high-aDC(A), pDC(B), CD8 $+$ T-cells(C), CD8 $+$ Tcm(D), CD8 $+$ Tem(E) xCell score in the GSE65904 cohort.

3.5 Prediction potential of the risk score in ICI therapy

The expression levels of PD-1, PD-L1, and CTLA-4 were significantly higher in the low-risk score group than the high-risk score group in the TCGA (Fig. 7A) and GSE65904 cohorts (Fig. 7B). Furthermore, patients with a low-risk score had a better survival outcome than those with a high-risk score in the TCGA (Fig. 7C, E, G) and GSE65904 cohorts (Fig. 7D, F, H).

Figure 7.

Impact of immune checkpoint gene expression and EMT-related genes signature on clinical outcome. Comparison of the expression pattern of immune checkpoint genes (PD1, PD-L1, and CTLA-4) between high-risk score group and low-risk score group in the TCGA cohort (A), GSE65904 cohort (B). Kaplan-Meier survival curves of OS among four patient groups stratified by the risk score and PD-1, PD-L1 and CTLA-4 in the TCGA cohort (C, E, G) and GSE65904 cohort (D, F, H).

3.6 Construction of nomograms

Age, race, T stage, N stage, and risk score were incorporated in the nomogram (Table 2). Moreover, a total point summarized the points of each variable, predicting the probability of overall survival at 1, 3, and 5 years (Fig. 8A). The area under curve (AUC) at 1, 3, and 5 years in the nomogram model showed higher accuracy (AUC $>$ 0.75) in overall survival (Fig. 8B). The calibration plots suggested that the nomogram performed well in comparison with the ideal model (Fig. 8C).

Table 2
Univariate and multivariate Cox analysis of TCGA cohort

Characteristics	Univariate Cox analysis			Multivariate Cox analysis
	HR	CI95%	$P$ value	HR	CI95%	$P$ value
Age	1.02	1.01–1.03	0.00008 ${}^{*}$	1.01	1–1.02	0.017 ${}^{*}$
Risk score	1.93	1.54–2.42	$P<$ 0.0001 ${}^{*}$	2	1.54–2.59	$P<$ 0.0001 ${}^{*}$
Sex	1.03	0.75–1.41	0.86601	–	–	–
Race	0.53	0.34–0.82	0.00448 ${}^{*}$	0.61	0.4–0.95	0.028 ${}^{*}$
T_stage	1.41	1.23–1.62	$P<$ 0.0001 ${}^{*}$	1.32	1.15–1.52	$P<$ 0.0001 ${}^{*}$
N_stage	1.47	1.27–1.7	$P<$ 0.0001 ${}^{*}$	1.53	1.31–1.79	$P<$ 0.0001 ${}^{*}$
M_stage	2.22	1.04–4.76	0.03957 ${}^{*}$	1.78	0.8–3.97	0.161
Morphology	0.89	0.64–1.23	0.48249	–	–	–
Metastasis	2.78	1.66–4.66	0.0001 ${}^{*}$	1.45	0.84–2.5	0.18
Neoadjuvant	0.71	0.36–1.39	0.31304	–	–	–

${}^{*}P<$ 0.05.

Figure 8.

Nomogram and correlation analysis. (A) Nomogram for predicting the probability of overall survival at 1, 3 and 5 years. (B) The 1, 3 and 5 years AUC of the nomogram model. (C) The 1, 3 and 5 years calibration plots of the nomogram.

Figure 9.

Immunohistochemical plots of the three hub genes associated with prognosis and statistical analysis of the positivity rate. (A, C, E) shows the protein expression of each gene in melanoma and in the paracancerous tissue. (B, D, F) shows the statistical analysis of the staining positivity rate for each gene in melanoma and in the paracancerous tissue. ${}^{*}$ Representative $P$ -value $<$ 0.05, ${}^{**}$ representative $P$ -value $<$ 0.01.

Figure 10.

Immunohistochemical plots of the four hub genes associated with prognosis and statistical analysis of the positivity rate. (A, C, E, G) shows the protein expression of each gene in melanoma and in the paracancerous tissue. (B, D, F, H) shows the statistical analysis of the staining positivity rate for each gene in melanoma and in the paracancerous tissue. ${}^{*}$ Representative $P$ -value $<$ 0.05, ${}^{**}$ representative $P$ -value $<$ 0.01.

3.7 Image analysis

All immunohistological images were observed under an inverted microscope and images were collected. Moreover, the staining differences between osteosarcoma and paraneoplastic tissue specimens were compared. The FBP1, STAT1, and PIM2 antibodies were significantly upregulated, whereas ESRP1, JAG2, NME1, and PSME3 were significantly downregulated in melanoma compared with that in paraneoplastic tissue. Representative pictures are showed in Figs 9 and 10. In addition, Fig. 9B, D, and F showed that the positive rate of immunohistochemical staining for FBP1, STAT1, and PIM2 was significantly higher in melanoma than that in paraneoplastic tissues. Fig. 10B, D, F, and H showed that the positive rate of immunohistochemical staining for ESRP1, JAG2, NME1, and PSME3 was significantly lower in melanoma than that in paraneoplastic tissues.

4. Discussion

The normal melanocytes can transition to metastatic melanomas because of EMT [14]. Moreover, EMT, which was reported to correlate with the tumor immune microenvironment, including the lymphoid and myeloid cells compartment, leads to antitumor and tumor-promoting effects by shifting the balance of the immune system [15]. EMT plays a crucial role in tumor progression because of the metastatic transition and immune system dysregulation.

We identified 28 EMT-related DEGs associated with immunity and metastasis to further reveal the relationship between EMT and tumor metastasis and TME. Finally, seven EMT-related genes out of the 28 DEGs were screened to establish the risk model, and the GSE65904 cohort was used for validation. The risk model showed a good prognosis and diagnosis predictive ability for melanoma. Moreover, the seven EMT-related genes that were incorporated in the risk model were associated with melanoma and its metastatic: Epithelial Splicing Regulatory Protein 1 (ESRP1), Signal Transducer and Activator of Transcription 1 (STAT1), Pim-2 Proto-Oncogene (PIM2), Jagged Canonical Notch Ligand 2 (JAG2), Fructose-Bisphosphatase 1 (FBP1), NM23 Nucleoside Diphosphate Kinase 1 (NME1), and Proteasome Activator Subunit 3 (PSME3). Moreover, JAG2 and NME1 can promote the expansion of cells with enhanced tumor and metastatic potential in the context of melanoma cell cultures [16, 17], whereas, FBP1, PIM2, and STAT1 act as a tumor suppressor in melanoma [18, 19, 20]. In addition, ESRP1 can work as informative biomarkers for immunotherapy because of the regulation ability of TIICs in melanoma [21]. Our results consolidated previous data that the EMT is strongly associated with melanoma and its metastatic potential.

Subsequently, we analyzed the relationship between EMT and TIICs in the TME. The results showed that three subtypes were downregulated, and 25 subtypes were upregulated in the low-risk group. The infiltration of immune cells, including aDC, pDC, CD8 $+$ T cells, and CD8 $+$ Tem, wherein the high infiltration of those prolongs the survival time of patients with melanoma, were significantly increased in the low-risk score group. Activated DC and pDC play a crucial role in antigen processing and presentation, thus initiating and shaping antitumor immune responses [22, 23]. CD8 $+$ T cells play critical roles in immunosurveillance against cancer [24]. Antigen uptake by DCs can be stored intracellularly and efficiently cross-presented to CD8 $+$ T cells, which play critical roles in tumor control, particularly in melanoma through direct cytotoxicity [25]. These results demonstrated that EMT in melanoma may induce infiltration of immune cells and that the interaction between them plays a key role in tumor progression.

PD-1, PD-L1, and CTLA-4, the common immune checkpoint molecules [26, 27, 28], were associated with the progression and prognosis of melanoma. ICI therapy has enormous potential in improving clinical outcomes for several tumors, including melanoma [29]. PD-1, PD-L1, and CTLA-4 are currently available biomarkers in clinical practice but not a sufficient independent predictor of ICI response [30]. Therefore, it is crucial to find predictive biomarkers that could identify patients who are most likely to benefit from ICI. In this study, the expression of immune checkpoint genes was much higher in the low-risk group. The survival analysis was performed for patients divided into four groups on the basis of the risk score and the immune checkpoint expression level. Moreover, we found that patients with melanoma at low risk and high immune checkpoint expression level have a better prognosis. These results indicated that ICI therapy might be more effective in patients in the low-risk group. Moreover, a previous study reported that ESRP1 can work as informative biomarkers for immunotherapy because its downregulation may promote the probability of responding to immune checkpoint blockade therapy to melanoma [31]. In addition, STAT1, an upstream gene in the pathway “PD-L1 expression and PD-1 checkpoint pathway in cancer”, was strongly associated with PD-L1 expression. These findings thus revealed that the risk model could work as a predictor for ICI response.

GSEA analysis was performed to further reveal the molecular mechanism between high- and low-risk groups. GSEA analysis showed that inflammatory response and immune-related pathways such as IL2/STAT5 signaling, interferon $\alpha$ / $\gamma$ response, cytokine-cytokine receptor interaction, and MAPK signaling pathway were significantly enriched in the low-risk score group. Moreover, antigen processing and presentation pathways, and apoptosis were strongly associated with a low-risk score. These findings suggested that the complex interplay among EMT, TIICs, and immune checkpoint gene in the TME has an effect on patient survival in melanoma, which is in line with the previous study in immunomodulatory interactions among EMT, immune cell infiltration, and checkpoint gene expression. Finally, the risk score was incorporated into the nomogram as an independent prognosis factor. Moreover, the nomogam can evaluate the prognosis of patients more accurately. However, more experiments are warranted to validate our findings.

5. Conclusion

A novel EMT-related risk score was developed and validated in an independent cohort. The individualized risk score could effectively conduct risk stratification, overall survival prediction, ICI therapy prediction, and immune microenvironment judgment for patients with melanoma, which would be conducive to patients’ precise treatment.

Declarations

Consent for publication

Not applicable.

Ethics approval and consent to participate

All subjects volunteered for the study and signed informed consent forms. In order to ensure confidentiality, the names of study participants were not included in the data. Information obtained from the data of the study participants is kept confidential. In addition, the Ethics Committee of The First Affiliated Hospital of Guangxi Medical University approved the study.

Data availability statement

The dataset from the TCGA database generated and/or analyzed during the current study are available in the TCGA dataset repository (https://tcga-data.nci.nih.gov/tcga/) and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). And all data and materials for this study shall be availed whenever requested by the editorial team, reviewers, and other users.

Competing interest

We declare that we do not have any competing interests.

Funding

This work was sponsored by the National Natural Science Foundation of China (81560359); National Natural Science Foundation of China (81860393). Funding bodies had no role in the study design, collection, analysis, and interpretation of the data or in writing the manuscript.

Author contributions

Zhan Xinli, Liang Tuo: Conceptualization, Methodology. Chen Jiarui, Xu GuoYong, Zhang Zide: Data curation, Investigation. Xue Jiang, Zeng Haopeng, Jiang Jie: Formal analysis, Software. Qin Zhaojie, Hao Li, Chen Tianyou, Ye Zhen, Nie Yunfeng: Visualization. Liang Tuo, Liu Chong: Writing- Reviewing and Editing.

Supplementary data

The supplementary files are available to download from http://dx.doi.org/10.3233/CBM-210329.

sj-xlsx-2-cbm-10.3233_CBM-210329.xlsx - Supplemental material

Supplemental material, sj-xlsx-2-cbm-10.3233_CBM-210329.xlsx

sj-xlsx-1-cbm-10.3233_CBM-210329.xlsx - Supplemental material

Supplemental material, sj-xlsx-1-cbm-10.3233_CBM-210329.xlsx

sj-pdf-1-cbm-10.3233_CBM-210329.pdf - Supplemental material

Supplemental material, sj-pdf-1-cbm-10.3233_CBM-210329.pdf

Footnotes

Acknowledgments

We are grateful to Dr. Xinli Zhan (Spine and Osteopathy Ward, The First Affiliated Hospital of Guangxi Medical University) for his kindly assistance in all stages of the present study.

References

Zhang

and Chou

, FGFR3 promotes the growth and malignancy of melanoma by influencing EMT and the phosphorylation of ERK, AKT, and EGFR, BMC Cancer 19 (2019), 963.

Ryu

S.H.

Heo

S.H.

Park

E.Y.

Choi

K.C.

Ryu

J.W.

Lee

S.H.

and Lee

S.W.

, Selumetinib Inhibits Melanoma Metastasis to Mouse Liver via Suppression of EMT-targeted Genes, Anticancer Res 37 (2017), 607–614.

Barbero

Castro

M.V.

Villanueva

M.B.

Quezada

M.J.

Fernández

N.B.

DeMorrow

and Lopez-Bergami

, An Autocrine Wnt5a Loop Promotes NF-κB Pathway Activation and Cytokine/Chemokine Secretion in Melanoma, Cells 8 (2019).

Song

Quach

Guo

Jang

G.B.

Maazi

Zhao

Sands

N.A.

Liu

G.K.

Peng

Yuan

Machida

Akbari

Hagiya

Yang

Punj

Tang

and Liang

, Transcriptional regulation of autophagy-lysosomal function in BRAF-driven melanoma progression and chemoresistance, Nat Commun 10 (2019), 1693.

Wang

Shen

Wei

and Li

, HIF-1α/VEGF signaling-mediated epithelial-mesenchymal transition and angiogenesis is critically involved in anti-metastasis effect of luteolin in melanoma cells, Phytother Res 33 (2019), 798–807.

Mao

X.Y.

Q.Q.

Gao

Y.F.

Zhou

H.H.

Liu

Z.Q.

and Jin

W.L.

, Gap junction as an intercellular glue: Emerging roles in cancer EMT and metastasis, Cancer Lett 381 (2016), 133–137.

Kuo

C.L.

Chou

H.Y.

Chiu

Y.C.

Cheng

A.N.

Fan

C.C.

Chang

Y.N.

Chen

C.H.

Jiang

S.S.

Chen

N.J.

and Lee

A.Y.

, Mitochondrial oxidative stress by Lon-PYCR1 maintains an immunosuppressive tumor microenvironment that promotes cancer progression and metastasis, Cancer Lett 474 (2020), 138–150.

Mahoney

K.M.

and Freeman

G.J.

, Learning from PD-1 Resistance: New combination strategies, Trends Mol Med 22 (2016), 448–451.

Taggart

Andreou

Scott

K.J.

Williams

Rippaus

Brownlie

R.J.

Ilett

E.J.

Salmond

R.J.

Melcher

and Lorger

, Anti-PD-1/anti-CTLA-4 efficacy in melanoma brain metastases depends on extracranial disease and augmentation of CD8(+) T cell trafficking, Proc Natl Acad Sci U S A 115 (2018), E1540–e1549.

10.

Wang

Saci

Szabo

P.M.

Chasalow

S.D.

Castillo-Martin

Domingo-Domenech

Siefker-Radtke

Sharma

Sfakianos

J.P.

Gong

Dominguez-Andres

W.K.

Mulholland

Azrilevich

Cordon-Cardo

Salmon

Bhardwaj

Zhu

and Galsky

M.D.

, EMT- and stroma-related gene expression and resistance to PD-1 blockade in urothelial cancer, Nat Commun 9 (2018), 3503.

11.

Zhao

Liu

Zheng

and Qu

, dbEMT 2.0: An updated database for epithelial-mesenchymal transition genes with experimentally verified information and precalculated regulation information for cancer metastasis, J Genet Genomics 46 (2019), 595–597.

12.

Lin

Huang

Luo

Zhu

Xiao

and Zhu

, Development of a prognostic index and screening of potential biomarkers based on immunogenomic landscape analysis of colorectal cancer, Aging (Albany NY) 12 (2020), 5832–5857.

13.

Jiang

Liu

Liang

Liao

Chen

Huang

Sun

Zhang

Wang

Chen

Yao

Chen

Guo

Liu

and Zhan

, TRIM68, PIKFYVE, and DYNLL2: The Possible Novel Autophagy- and Immunity-Associated Gene Biomarkers for Osteosarcoma Prognosis, Front Oncol 11 (2021), 643104.

14.

Soundararajan

Fradette

J.J.

Konen

J.M.

Moulder

Zhang

Gibbons

D.L.

Varadarajan

, Wistuba, II Tripathy

Bernatchez

Byers

L.A.

Chang

J.T.

Contreras

Lim

Parra

E.R.

Roarty

E.B.

Wang

Yang

Barton

Rosen

J.M.

and Mani

S.A.

, Targeting the interplay between epithelial-to-mesenchymal-transition and the immune system for effective immunotherapy, Cancers (Basel) 11 (2019).

15.

Kudo-Saito

Shirako

Takeuchi

and Kawakami

, Cancer metastasis is accelerated through immunosuppression during Snail-induced EMT of cancer cells, Cancer Cell 15 (2009), 195–206.

16.

Asnaghi

Handa

J.T.

Merbs

S.L.

Harbour

J.W.

and Eberhart

C.G.

, A role for Jag2 in promoting uveal melanoma dissemination and growth, Invest Ophthalmol Vis Sci 54 (2013), 295–306.

17.

Wang

Leonard

M.K.

Snyder

D.E.

Fisher

M.L.

Eckert

R.L.

and Kaetzel

D.M.

, NME1 drives expansion of melanoma cells with enhanced tumor growth and metastatic properties, Mol Cancer Res 17 (2019), 1665–1674.

18.

Gutteridge

R.E.

Singh

C.K.

Ndiaye

M.A.

and Ahmad

, Targeted knockdown of polo-like kinase 1 alters metabolic regulation in melanoma, Cancer Lett 394 (2017), 13–21.

19.

Chen

Crabill

G.A.

Pritchard

T.S.

McMiller

T.L.

Wei

Pardoll

D.M.

Pan

and Topalian

S.L.

, Mechanisms regulating PD-L1 expression on tumor and immune cells, J Immunother Cancer 7 (2019), 305.

20.

Shannan

Watters

Chen

Mollin

Dörr

Meggers

Gimotty

P.A.

Perego

Benci

Krepler

Brafford

Zhang

Wei

Zhang

Liu

Yin

Nathanson

K.L.

Herlyn

and Vultur

, PIM kinases as therapeutic targets against advanced melanoma, Oncotarget 7 (2016), 54897–54912.

21.

Yao

Caballero

O.L.

Huang

Lin

Rimoldi

Behren

Cebon

J.S.

Hung

M.C.

Weinstein

J.N.

Strausberg

R.L.

and Zhao

, Altered Expression and Splicing of ESRP1 in Malignant Melanoma Correlates with Epithelial-Mesenchymal Status and Tumor-Associated Immune Cytolytic Activity, Cancer Immunol Res 4 (2016), 552–561.

22.

Aspord

Leccia

M.T.

Charles

and Plumas

, Plasmacytoid dendritic cells support melanoma progression by promoting Th2 and regulatory immunity through OX40L and ICOSL, Cancer Immunol Res 1 (2013), 402–415.

23.

Aspord

Tramcourt

Leloup

Molens

J.P.

Leccia

M.T.

Charles

and Plumas

, Imiquimod inhibits melanoma development by promoting pDC cytotoxic functions and impeding tumor vascularization, J Invest Dermatol 134 (2014), 2551–2561.

24.

Blanch-Lombarte

Gálvez

Revollo

Jiménez-Moyano

Llibre

J.M.

Manzano

J.L.

Boada

Dalmau

E.S.

Clotet

Prado

J.G.

and Martinez-Picado

, Enhancement of Antiviral CD8(+) T-Cell Responses and Complete Remission of Metastatic Melanoma in an HIV-1-Infected Subject Treated with Pembrolizumab, J Clin Med 8 (2019).

25.

Edwards

Wilmott

J.S.

Madore

Gide

T.N.

Quek

Tasker

Ferguson

Chen

Hewavisenti

Hersey

Gebhardt

Weninger

Britton

W.J.

Saw

R.P.M.

Thompson

J.F.

Menzies

A.M.

Long

G.V.

Scolyer

R.A.

and Palendira

, CD103(+) Tumor-Resident CD8(+) T Cells Are Associated with Improved Survival in Immunotherapy-Naïve Melanoma Patients and Expand Significantly During Anti-PD-1 Treatment, Clin Cancer Res 24 (2018), 3036–3045.

26.

Luo

Firat

Gaedicke

Guffart

Watanabe

and Niedermann

, Cisplatin Facilitates Radiation-Induced Abscopal Effects in Conjunction with PD-1 Checkpoint Blockade Through CXCR3/CXCL10-Mediated T-cell Recruitment, Clin Cancer Res 25 (2019), 7243–7255.

27.

Donia

Kjeldsen

J.W.

Andersen

Westergaard

M.C.W.

Bianchi

Legut

Attaf

Szomolay

Ott

Dolton

Lyngaa

Hadrup

S.R.

Sewell

A.K.

and Svane

I.M.

, PD-1(+) Polyfunctional T Cells Dominate the Periphery after Tumor-Infiltrating Lymphocyte Therapy for Cancer, Clin Cancer Res 23 (2017), 5779–5788.

28.

Wei

S.C.

Levine

J.H.

Cogdill

A.P.

Zhao

Anang

N.A.S.

Andrews

M.C.

Sharma

Wang

Wargo

J.A.

Pe’er

and Allison

J.P.

, Distinct Cellular Mechanisms Underlie Anti-CTLA-4 and Anti-PD-1 Checkpoint Blockade, Cell 170 (2017), 1120-1133e17.

29.

Loo

and Daud

, Emerging biomarkers as predictors to anti-PD1/PD-L1 therapies in advanced melanoma, Immunotherapy 8 (2016), 775–784.

30.

Ahern

Harjunpää

Barkauskas

Allen

Takeda

Yagita

Wyld

Dougall

W.C.

Teng

M.W.L.

and Smyth

M.J.

, Co-administration of RANKL and CTLA4 antibodies enhances lymphocyte-mediated antitumor immunity in mice, Clin Cancer Res 23 (2017), 5789–5801.

31.

Wang

Kou

Sun

and Xu

, Mining Database for the Clinical Significance and Prognostic Value of ESRP1 in Cutaneous Malignant Melanoma, Biomed Res Int 2020 (2020), 4985014.

Epithelial-mesenchymal transition interaction with CD8 + T cell,dendritic cell and immune checkpoints in the development of melanoma

Abstract

BACKGROUND:

OBJECTIVE:

METHODS:

RESULTS:

CONCLUSION:

Keywords

Abbreviations

1. Introduction

2. Patients and methods

2.1 Microarray data acquisition

2.2 Identification of DEGs and PPI network construction

2.3 Establishment of a prognostic gene signature for melanoma

2.4 Gene set enrichment analysis (GSEA)

2.5 xCell score estimation

2.6 Relationships between risk score, immune cells and immune checkpoints

2.7 Construction of nomograms

2.8 Immunohistochemistry

Table 1 Univariate Cox regression analyses to the DEGs

3.1 EMT-related DEGs were mainly involved in inflammation response

3.3 The risk score predicts the involvement of immune pathways

3.4 TIICs profile in different risk groups

Table 2 Univariate and multivariate Cox analysis of TCGA cohort

4. Discussion

5. Conclusion

Declarations

Consent for publication

Ethics approval and consent to participate

Data availability statement

Competing interest

Funding

Author contributions

Supplementary data

sj-xlsx-2-cbm-10.3233_CBM-210329.xlsx - Supplemental material

sj-xlsx-1-cbm-10.3233_CBM-210329.xlsx - Supplemental material

sj-pdf-1-cbm-10.3233_CBM-210329.pdf - Supplemental material

Footnotes

Acknowledgments

References

Table 1
Univariate Cox regression analyses to the DEGs

Table 2
Univariate and multivariate Cox analysis of TCGA cohort