CD247 Functions as a Prognostic Biomarker for Cutaneous Malignant Melanoma Based on the Analysis of Tumor-immune Microenvironment

doi:10.21203/rs.3.rs-54682/v1

Download PDF

Research article

CD247 Functions as a Prognostic Biomarker for Cutaneous Malignant Melanoma Based on the Analysis of Tumor-immune Microenvironment

https://doi.org/10.21203/rs.3.rs-54682/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Cutaneous malignant melanoma (CMM) is among the most lethal cancers. The tumour microenvironment (TME) is closely linked with tumorigenesis, metastasis, and prognosis.

Methods: We employed the ESTIMATE algorithm to calculate immune and stromal scores of malignant melanoma tissues from the Cancer Genome Atlas dataset, respectively, and determine core prognosis gene signature examined by COX proportional hazards model. Functional enrichment annotation, the Kyoto Encyclopedia of Genes and Genomes pathways, the Protein-Protein Interaction network, Weighted Gene Co-expression Network Analysis and overall survival analysis were used to formulate potential function of these genes that involved in immune-linked biological processes. CIBERSORT algorithm was used to estimate the abundances of immune cell types in CMM samples. Finally, the OncoLnc platform, the Gene Expression Profiling Interactive Analysis resources, and the Human Protein Atlas database were applied to validate our results.

Results：908 differential expressed genes and ten hub genes were screened, and GO annotation indicated that immune response and inflammatory response were firmly involved in CMM tumorigenesis and progression. CD247, identified as the most significant prognostic biomarker, highly expressed in tumor samples and possessed a better prognosis than low expressed samples. The correlation analysis of immune cells infiltration unveiled that CD8⁺ T cell and Macrophages were intense significant to CMM patients' prognosis. Survival analysis suggested that ten hub genes and infiltrated immune cells are linked to the prognosis of CMM. ConnectiveMap analysis strongly indicated that L-securinine may be a promising candidate medicine for CMM patients.

Conclusions: we deeply analyzed the immune-linked genes with the tumour microenvironment, and labeled CD247 as the most intriguing prognostic biomarker for CMM, which may bring better clinical outcomes for CMM patients.

Cancer Biology

Oncology

tumour-microenvironment

bioinformatics

cutaneous melanoma

prognostic biomarker

Cutaneous malignant melanoma (CMM), one of the most aggressive cancer with a reduced survival rate, features the primary cause of skin cancer–related death, and is notorious for its resist-therapy[1]. In the past several years, immunotherapeutic strategies have made significant progress and did improve patient survival time. However, owing to CMM's easy-to-recurrence and dreadful metastasis, amounts of patients adapting such treatments have no obvious durable response[2, 3]. Besides, some adverse effects have emerged, such as autoimmunity, since the complexity of interaction between tumour cells and the tumour microenvironment (TME). TME comprises immune cells, endothelial cells, mesenchymal cells, inflammatory factors, and extracellular matrix (ECM) molecules[4]. The cells and particles in the TME are in a dynamic process and supported to be highly linked with various tumour behaviours, including proliferation, metastasis, and immune escape. Immune cells and stromal cells are the two main types of non-tumour components and considered to be of significant importance in the diagnosis and prognosis of tumours [5]. An increasing amount of evidence has elucidated the clinic pathological significance of TME in the prediction of treatment effects[6]. Therefore, understanding the molecular composition and function of TME has more effective regulation of cancer progression and immune response in CMM.

Extensive effort has been made on the mechanism of the occurrence and progression of CMM. However, the concrete pathogenesis of CMM still needs to be further formulated. Therefore, it is an urgent task to explore freshly molecular markers that represent a diagnostic and prognostic value for CMM. In this study, we comprehensively investigated the TME and immune cells associated with CMM to determine immune-related prognostic molecules for CMM. Finally, the HPA database, the OncoLnc platform, and the GEPIA tool were applied to validate the study outcomes. The whole experiment procedure was shown in figure 1.

2.1 | Raw data collection

The RNA-seq expression profiling for CMM patients was acquired from the TCGA database (https://tcga‐data.nci.nih.gov/tcga/gdc), including 32 normal samples and 436 tumour samples[7]. The workflow type was selected HTSeq-FPKM for further investigation, which including 472 files and 468 cases. Clinical data involving each patients’ age, TNM staging, gender, tumour grade, and survival information were also downloaded from the TCGA website. The gene expression quantification was measured using the Illumina HiSeq2000 RNA Sequencing platform. Then the ESTIMATE algorithm (http://r-forge.r-project.org) was applied to calculate each samples’ immune and stromal scores, respectively, which was performed by R software (4.0.2) with the help of relevant R packages: “estimate,” “limma” and “utils.” To validate our results[8], we explored the OncoLnc website (http://www.oncolnc.org/), GEPIA platform (http://gepia.cancer-pku.cn/), and the HPA (https://www.proteinatlas.org/) to verify the screened prognostic-genes.

2.2 | Differential analysis of expressed genes

Initially, all samples were separated into high and low immune‐score groups and high and low stromal‐score groups based on the outcomes of the ESTIMATE analysis. Then, the R package “Limma” was used to identify the differential expression genes (DEGs) between high and low immune‐score groups, as well as high and low stromal‐score groups[9]. We set the screening conditions as the DEGs were: log₂ |fold change| > 1.5 and adjust p-values <0.05. Heatmaps of DEGs were drawn via the R package “pheatmap.”

2.3 | WGCNA construction

WGCNA was conducted by adopting the R packages of “wgcna,” matrixStats,” dynamic tree cut,” fast cluster,” Hmisc” in R software[10]. Pearson’s correlation matrices were prepared for all pair-wise genes. Then, a power function amn =ⅠcmnⅠβ (cmn = Pearson’s correlation between gene m and n; amn = adjacency between gene m and n) was selected to erect weighted adjacency matrix. We transformed the adjacency to the topological overlap matrix (TOM), which tests the network connectivity of genes. Then, we developed average linkage hierarchical clustering in line with the TOM-based dissimilarity measure with a minimum size of 50 for the genes dendrogram. Therefore, genes that possessed similar expression profiles might be classified into the same module, which defined as groups of genes with a high degree of correlation. Then a single‐column vector called the module eigengenes were built to elucidate the potential relationship between the gene modules and clinicopathological parameters. The module eigengenes were representative of the gene expression profiles in each module, which was produced by preserving the first principal component of a designated module.

Moreover, the module eigengenes represent a summary metric for the co-expression genes network for each module eigengene contains most of the variance in the raw data. The consistency between the module eigengenes expression and the expression of genes was labelled as the module membership. We calculated the dissimilarity of module eigengenes, selecting a cut line for module dendrogram to merge them. Functional analysis was performed based on the modules, which have a notable influence on CMM.

2.4 | Functional Enrichment annotation and PPI network

Functional enrichment analysis was performed by using the R packages: “clusterProfiler,” “enrichplot,” “org.Hs.eg.db,” “ ggplot2” in R software to annotation, visualization and integrated outcomes of the DEGs, the corresponding biological processes (BP), cell components (CC), and molecular functions (MF) were determined via Gene Ontology (GO)[11], and the most enriched signalling pathways were identified through the Kyoto Encyclopedia of Genes and Genomes (KEGG)[12]. A false discovery rate < 0.05 was considered as the cut‐off. The protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/cgi/input.pl), then, the top ten hub genes were identified using the CytoHubba, plug-in of Cytoscape (3.8.0) software. Module analysis for the detection of interaction networks was performed using the MCODE plug-in of Cytoscape platform[13].

2.5 | Clinical correlation analysis of CD247

For determining the core prognostic gene, the R package “VennDiagram” was utilized. We combined the DEGs obtained from PPI network and prognostic relevant genes acquired from the analysis of the COX regression model. Then, CD247 represented the most significant prognostic gene. Next, a series of R packages: “limma,” “beeswarm,” “ggpubr,” “survival,” and “survminer” were employed in R software to investigate the clinical correlation of CD247. More critical, the GEPIA tool (http://gepia.cancer-pku.cn/) was applied to verify the prognostic value of CD247 in other CMM cases.

2.6 | Immune cells infiltration analysis

To better reveal that the underlying relationship between infiltrating immune cells, including CD8⁺ T cell, CD4⁺ T cell, B cells, eosinophils, neutrophils, macrophages, and dendritic cell in CMM, and prognosis in CMM patients, R software was applied to analyze the clinical correlation of infiltrating immune cells with the help of R packages: “limma,” “preprocesscore,” “corrplot” and “vioplot.” CIBERSORT algorithm, an analytical tool, to provide an estimation of the abundances of immune cell types in a mixed cell population performed by disease gene expression matrix[14].

2.7 | ConnectiveMap (CMap) and overall survival analysis

CMap (https://portals.broadinstitute.org/cmap/) is an open database that association with disease, genes, and drugs or small compounds based on the gene expression profile. Therefore, CMap analysis was used to identify potential drugs or small compounds to mitigate cutaneous malignant melanoma. We regarded Mean < −0.4 and p-value < 0.05 as the screening standard. The survival analysis was shown using the Kaplan–Meier curve with a log-rank test, which is conducted by the R packages “survival” and “survminer.” The survival curve illustrates the association between prognostic-genes and overall survival time.

2.8 | Statistical methods

All data were demonstrated as mean ± SD. One‐way analysis of variance was used to compare the stromal and immune scores in assigned groups by GraphPad Prism 8.0 software, and a two-tailed p < 0.05 was regarded as significant.

3.1 | Stromal and immune scores are highly linked with the TNM staging system and prognosis

The RNA-seq profiling and clinicopathologic information about 468 patients with CMM were downloaded from the TCGA database. And then, we processed and analyzed the samples based on the workflow (see Fig. 1). Among all samples, the age of individuals was 62\(\pm\)8. 95.1% were white (n=445), 2.6% Asian (n=12), 2.3% were not reported (n = 10) and others were black or African–American. The proportion of patients with T1–T2 and T3–T4 was 37.6% (n = 176) and 62.4% (n = 292), respectively. Patients with N0–N1 and N2–N3 stand at 66.2% (n = 310) and 33.8% (n = 158), respectively. Patients with M0 and M1 reached 87.2% (n = 408) and 12.8% (n = 60), respectively (see additional file 3). Next, the ESTIMATE algorithm was applied to calculate stromal scores and immune scores of all samples, and the score was ranged from − 1806.847 to 1891.958 and − 1481.104 to 3769.121, respectively (see additional file 1). As for the N staging and M staging of CMM patients, both immune score and stomal score reveal no significant association between N/M staging and immune/stromal scores (Fig. 2e, 2f, and additional file 2). In terms of T staging, the order for stromal median score is: T4 < T3 < T2 < T1, and order for immune median score represent: T4 < T2 < T3 <T1. There is a significant difference between them (Fig. 2d). To unveil the underlying relationship between the stromal/immune/estimate scores and the overall survival (OS) of the CMM patients, we divided samples into high and low score groups according to the high/low stromal/ immune/estimate scores. Then, the Kaplan–Meier survival curve showed that the high score group of the immune scores has a higher survival rate than the low score group (Fig. 2a, p<0.001). Similar outcomes were observed in the high and low score groups of the immune/estimate scores (Fig. 2b, p = 0.076, Fig. 2c, p < 0.001).

3.2 | Analysis of Differential expression genes (DEGs) from CMM cases based on stromal and immune scores

In our current study, differential analysis of all RNA-seq data from 468 CMM cases was performed to understand better the relationship between the stromal/immune scores and the gene expression profile of the samples (see additional file 4 and 5). Then, according to the distribution of the stromal/immune scores, a heatmap was drawn to reveal DEGs profiles of the samples, in which 1352 upregulated genes and 36 down-regulated genes were acquired based on the difference in stromal scores (log₂FC > |1.5|, p-adj < 0.05). Similarly, 1042 upregulated genes and 94 down-regulated genes were obtained based on the differential analysis of immune scores (log₂FC > |1.5|, p-adj < 0.05) (Fig. 3a). Notably, immune-related genes were explored via the intersection genes that were combined the upregulated or downregulated in both immune and stromal groups (Fig. 3b. 900 upregulated genes and eight downregulated genes, see additional file 6).

3.3 | Weighted gene co-expression network analysis

The R package “WGCNA” was loaded in R software and set the power β to a soft-thresholding parameter β = 6 (scale-free R² = 0.87). MEDissThres was defined as 0.25 to acquisition similar modules, and then 12 modules were generated. The black module included 221 genes, and the blue module included 1935 genes. While the brown module contained 2365 genes, and the green module contained 675 genes. Grey module represented 721 genes. The red module included 239 genes. The turquoise module included 1989 genes. Greenyellow module included 138 genes. The magenta module included 160 genes. The pink module included 220 genes. The purple module included 148 genes. Genes in the yellow module did not belong to any of the functional modules, and we neglected it. Besides, an intramodular analysis of Gene significance (GM) and module membership (MM) of the genes in the 11 modules involved was performed. Since GS and MM illustrated meaningful correlation, the result implied that the brown module was strongly related to CMM among the 11 modules. Moreover, genes belong to module brown had the highest positive correlation with CMM (see additional file 7). In contrast, genes in the red module had the highest negative ties with CMM (see additional file 8).

3.4 | Functional Enrichment annotation and PPI network analysis of DEGs

We predicted the function of the intersection 908 differential expression genes (900 upregulated and eight down-regulated). Go annotation (GO), and KEGG pathway analysis was performed by R software, including biological process (BP), molecule function (MF), cellular component (CC), and enriched signalling pathways[15]. Sorting by adjusting the p-value, we select the top 10 terms of each section (see Table 1). The results unveiled that GO functions are mainly enriched in inflammatory and immune responses, T cell activation, regulation of lymphocyte activation, and regulation of T cell activation (Fig. 5a, 5c). In comparison, KEGG pathways are mainly enriched in Cytokine-cytokine receptor interaction and Hematopoietic cell lineage (Fig. 5b, Table 2). To analyze the association between genes with prognostic value, the STRING network tool was used to construct a PPI network of these genes (Fig. 6a, 6b). Hub genes analysis of the PPI network was performed via the CytoHubba, a plug-in in Cytoscape, the core gene scores were calculated using the Maximal Clique Centrality (MCC) method[16]. Then, the top 10 key genes were shown as follows: CCL5, CXCL10, CD74, CXCL9, CD247, IL10, CXCL11, CXCR3, CD3E and SYK (Fig. 6c, Table 2). The module contains 26 nodes and 67 edges. Cytoscape’s MCODE was used to conduct a comprehensive modular analysis of the differential genes (Fig. 6d). Three modules based on scores were identified through MCODE analysis, and expanded nodes and edges were not displayed. (screen criteria: Degree cut-off: 2, Node score cut-off: 0.4, K-core: 2). The module one contains 8 points and 16 edges with a 4571 score, which including CCL5, CCL3, CCR5, CCL4, CXCR3, CXCL9, CXCL10, CXCL11. In module two, with five nodes and ten edges marked 4500 scores, which containing CYBB, NCF1, RAC2, NCF2, NCF4. In module three, with six nodes and nine edges marked 4000 scores, which possessing CD3E, LCK, CD3G, ZAP70, CD3D, and CD28 (Fig. 6). Next, the above-talked genes identified by PPI key modules functional analysis were investigated using GO annotation and KEGG pathway in R software and visualized by ClueGo, a plug-in in Cytoscape (Kappa score: 0.4, Go tree interval: min-level:4, max-level:8, and p-value < 0.05). The outcomes demonstrated that these genes were largely involved in the immune and inflammatory response, regulation of lymphocyte activation, leukocyte proliferation and play a critical role in Cytokine-cytokine receptor interaction, Chemokine signalling pathway and cell adhesion molecules (CAMs) (Fig. 6e).

Table 1

List of GO annotation and KEGG enrichment pathway of DEGs.
Ontology	ID	Description	p-adjust	q-value	Count
BP	GO:0042110	T cell activation	2.68E-80	2.01E-80	145
BP	GO:0051249	regulation of lymphocyte activation	1.82E-68	1.36E-68	136
BP	GO:0007159	leukocyte cell-cell adhesion	9.12E-68	6.83E-68	115
BP	GO:0050863	regulation of T cell activation	3.78E-66	2.83E-66	110
BP	GO:0070661	leukocyte proliferation	4.06E-61	3.04E-61	103
MF	GO:0140375	immune receptor activity	2.80E-34	2.41E-34	53
MF	GO:0004896	cytokine receptor activity	5.89E-17	5.07E-17	32
MF	GO:0005126	cytokine receptor binding	9.23E-17	7.94E-17	54
MF	GO:0005125	cytokine activity	9.23E-17	7.94E-17	47
MF	GO:0008009	chemokine activity	6.36E-15	5.47E-15	22
CC	GO:0009897	external side of plasma membrane	4.21E-45	3.31E-45	99
CC	GO:0030667	secretory granule membrane	2.64E-18	2.08E-18	56
CC	GO:0042613	MHC class II protein complex	1.25E-15	9.82E-16	14
CC	GO:0070820	tertiary granule	4.15E-13	3.27E-13	35
CC	GO:0042611	MHC protein complex	6.32E-13	4.98E-13	15
KEGG	hsa04060	Cytokine-cytokine receptor interaction	4.64E-33	3.76E-33	84
KEGG	hsa04640	Hematopoietic cell lineage	1.29E-31	1.04E-31	49
KEGG	hsa04061	Viral protein interaction with cytokine and cytokine receptor	1.59E-31	1.29E-31	49
KEGG	hsa04672	Intestinal immune network for IgA production	6.77E-23	5.49E-23	30
KEGG	hsa04514	Cell adhesion molecules (CAMs)	2.42E-21	1.97E-21	48
BP: biological process, MF: molecular function, CC: cell component
DEGs: differential expression genes.

Table 2

Characteristics and function role of ten hub genes.
No.	Gene	Full Name	Function
1	CCL5	C-C chemokine ligand 5	CCL5 promotes tumour progression and recruitment of immune infiltration leukocytes in several cancers and plays a vital role in inflammatory diseases.
2	CXCL10	CXC motif chemokine ligand 10	CXCL10 involved in tumour cell invasion and migration by activating the TLR signalling pathway with inhibiting angiogenesis and influence cancer regression
3	CD74	cluster of differentiation 74	CD74, a cell membrane receptor for macrophage migration inhibitory factor (MIF), regulates T-cell and B-cell developments, dendritic cell (DC) motility, macrophage inflammation, and thymic selection.
4	CXCL9	CXC motif chemokine ligand 9	CXCL9 binding to a subset of G protein-coupled receptors and play a crucial role to induce chemotaxis, promote differentiation and multiplication of leukocytes, as well as cause tissue extravasation.
5	CD247	Cluster of differentiation 247	CD247 activated the T cell receptors (TCR) signalling cascade and promote assembling of the TCR/CD3 complex on the surface of T lymphocytes.
6	IL10	Interleukin-10	IL10 promotes CD8⁺T cell cytolytic function, inhibits antigen presentation, and production of proinflammatory cytokine.
7	CXCL11	CXC motif chemokine ligand 11	CXCL9, CXCL10, and CXCL11 are core immune chemo-attractants during inflammatory responses, and CXCL11 attenuates tumour angiogenesis.
8	CXCR3	CXC motif chemokine receptor 3	CXCR3 is an interferon-inducible chemokine receptor and mediates functions of CXC motif chemokine ligands(CXCL9, CXCL10, CXCL11).
9	CD3E	Cluster of differentiation 3ε	CD3E is one part of TCR, involved in antigen recognition.
10	SYK	Spleen tyrosine kinase	SYK, mainly expressed in hematopoietic cells, activates the B-cell receptor signalling pathway and mobilizes calcium ion channel.

3.5 | CD247 was determined as the most valuable prognostic gene.

For a thoroughly investigated the relationship between identified prognostic-involved genes and overall survival in CMM patients, a univariate Cox model was employed to select prognostic-related genes highly correlated with overall survival (see Table 3, Fig. 7d). Next, we combined the prognostic-related genes with the immune-related DEGs from the PPI network, and the result showed that CD247 was the intersection gene. Therefore, CD247 was regarded as the critical prognostic-related gene for further study. To begin with, we compared the expression of CD247 in normal cases (n = 32) vs tumour cases (n = 436), and the result showed that CD247 highly expressed in tumour samples (p = 0.011) and high expression of CD247 cases have a higher survival rate than low group (p < 0.001). The expression of CD247 in both female patients and male patients have no significant difference (p = 0.12) (Fig. 7a). Then, the analysis of TNM staging was performed to support the prognostic value of CD247 (Fig. 7b), from the perspective of tumour infiltration depth (T staging), the expression of CD247 from T1 vs T2 and T1 vs T3 to T1 vs T4 demonstrated a sharp significant difference. At the same time, there was no meaningful relationship between CD247 and N/M staging. (Fig. 7b). Also, we mined the GEPIA online platform and collected 1019 CMM cases, which covered 461 tumour samples (n = 461) and 558 normal samples (n = 558). Then, we inspected the expression of CD247 in each sample between tumour and normal tissues as well as performed survival curve to convince prognostic characterize of CD247. As is showed in Fig. 8c, CD247 is higher expression in tumour samples compared with normal samples (p < 0.05). Survival analysis revealed that the CD247 higher expression group possesses a better prognosis than the lower group (Fig. 7c).

Table 3

Prognostic-linked genes produced by COX proportional hazards model.
gene	KM	HR	HR.95L	HR.95H	p-value
APOL3	1.39E-07	0.913736031	0.882228585	0.946368717	4.68E-07
P2RY6	4.61E-05	0.599792304	0.487871077	0.737389086	1.23E-06
LILRB1	1.37E-06	0.780246359	0.703119834	0.865833037	2.97E-06
CD72	3.00E-06	0.829776761	0.766191293	0.898639125	4.49E-06
HK3	8.76E-07	0.819926255	0.752899264	0.892920335	5.05E-06
LAG3	1.73E-06	0.936318983	0.910037577	0.963359383	5.90E-06
IL15RA	1.54E-08	0.767115381	0.683103301	0.861459763	7.47E-06
IL10RA	9.86E-07	0.913321751	0.877534643	0.950568309	8.76E-06
C19orf38	1.61E-05	0.749694891	0.655296348	0.85769199	2.72E-05
CCL5	9.47E-08	0.991825762	0.988013085	0.995653152	2.96E-05
SOCS1	0.000165132	0.872088961	0.817763228	0.930023666	3.04E-05
FYB1	9.78E-07	0.937153531	0.90866344	0.966536895	3.78E-05
RGL4	7.59E-05	0.119698237	0.043415569	0.330012213	4.09E-05
CD247	4.13E-06	0.216395217	0.103892046	0.450726423	4.34E-05
LRRC25	5.22E-06	0.891676204	0.842416454	0.943816385	7.68E-05
RTP5	8.89E-08	0.057189895	0.013748227	0.237898608	8.34E-05
TBC1D10C	2.24E-05	0.862279301	0.799567112	0.929910175	0.000119977
ZBP1	6.16E-06	0.643455126	0.513394977	0.806463867	0.000129735
FBP1	5.13E-06	0.945042313	0.917860198	0.973029417	0.000146982
FLT3LG	0.000210684	0.473961697	0.321736036	0.698211158	0.000158442
CASS4	9.87E-06	0.403369568	0.245805838	0.661933049	0.000327394
CEACAM4	0.003479028	0.240665423	0.109753769	0.527725349	0.000377253
IRF5	0.002693029	0.850900132	0.775098771	0.934114544	0.000694697
SIGLEC5	0.002679193	0.09131943	0.022151131	0.376470089	0.00092726
DBH	0.000398379	0.192496102	0.072026844	0.514457489	0.001019398
ADAMDEC1	1.56E-05	0.937700212	0.902345163	0.974440517	0.001036749
C4B	0.004414878	0.688641099	0.540152816	0.877948886	0.002608862
C11orf21	0.015153385	0.539364469	0.35944481	0.809342693	0.002867803
OTULINL	9.11E-08	0.863808702	0.783449693	0.952410193	0.00329607
ICAM2	0.006961249	0.826749292	0.72665423	0.940632235	0.003858699
XIRP1	2.08E-05	0.535727612	0.339677767	0.844930408	0.007257321
SERPING1	2.90E-05	0.996898696	0.994615389	0.999187243	0.007932087
ASGR2	0.000704977	0.38601548	0.187731321	0.793729837	0.009651716
SSTR3	1.25E-05	0.189201927	0.050838537	0.704138455	0.013023556
CEBPA	0.000489433	0.926792035	0.872650051	0.98429316	0.013307022
DPEP2	0.001484138	0.757579093	0.60493555	0.94873922	0.015591784
STAB1	0.013789578	0.956505019	0.920736465	0.993663102	0.02220263
SELL	5.51E-05	0.983351676	0.968610095	0.998317614	0.029372381
KM: Kaplan–Meier. HR: hazard ratio.

3.6 | The relationship between immune cells infiltration abundance and prognosis of CMM

The genes expression matrix and clinical information were applied to assess the abundance of immune infiltrating cells in each case based on CIBERSORT algorithm and to unveil the relationship between immune cells infiltration and CMM prognosis (Fig. 8a, 8c). The vioplot showed the distribution of immune cells. Among them, the highest proportion of immune cells was Macrophage M0, followed by CD8⁺ T cells (Fig. 8b). Then, the prognostic correlation analysis was performed for immune cell infiltration, and CD247 was drawn (Fig. 8d). The results showed that all of which were functioned as a critical role in CMM microenvironment and equipped with prognostic values for CMM patients’ survival time, among them, MacrophagesM0, MacrophagesM2, rested CD4 memory T cell and Eosinophils had a negative correlation with CD247 while CD8⁺ T cell, MacrophagesM1, activated CD4 memory T cell, regulatory T cell, helper follicular T cell, and activated NK cell has a positive correlation with CD247 (see Table 5).

3.7 | The survival curve of prognostic-related hub genes and ConnectiveMap (CMap) analysis

the Kaplan–Meier curve was drawn to visualize survival analysis, which is performed using R packages: “survival” and “survminer” in R software. A log-rank test was to verify the association between the ten key genes and the overall survival of CMM cases. As is demonstrated in Fig. 9. Exploring the CMap platform, we identified drugs or small compounds that may help ameliorate CMM. Mean < − 0.4 and p < 0.05 as the screening condition (see Table 4).

Table 4

Drugs or small molecules explored by CMap analysis with mean < -0.4 and P < 0.05.
cmap name	mean	enrichment	p	per cent non-null
STOCK1N-35696	-0.765	-0.954	0.00451	100
securinine	-0.641	-0.728	0.01122	75
harmol	-0.57	-0.712	0.0141	75
harpagoside	-0.549	-0.726	0.01154	75
dihydroergocristine	-0.547	-0.801	0.00308	75
atractyloside	-0.466	-0.692	0.00619	60
gentamicin	-0.46	-0.833	0.00143	50
dihydrostreptomycin	-0.41	-0.663	0.01029	60

Table 5

Response p-value of immune cells.
Cell	p-value
T cells CD8	1.47E-51
Macrophages M0	7.01E-19
T cells CD4 memory activated	1.30E-16
T cells CD4 memory resting	1.27E-09
T cells follicular helper	6.07E-08
Macrophages M1	4.95E-07
Plasma cells	2.14E-05
T cells gamma delta	0.0001186
Dendritic cells resting	0.0001815
Eosinophils	0.0002302
Macrophages M2	0.0002532
NK cells resting	0.0008666
NK cells activated	0.0069879
T cells regulatory (Tregs)	0.012249

3.8 | Mining the OncoLnc database and HPA platform to validate obtained prognostic genes

The OncoLnc platform (http://www.oncolnc.org/) was used to verify whether the ten prognostic related genes in the TCGA database also crucial for other CMM cases. We investigated 458 samples (n = 458) and applied the Kaplan–Meier curve to analyze survival time for CMM patients (Fig. 10). The results revealed that the higher expression of these genes may be potential indicator genes for a good prognosis and may offer new insight into the therapies for CMM. To further confirm the reliability of the prognostic-response genes, IHC evidence from the HPA database was employed to display the expression of these genes (CD74, CD247, IL10, CXCL11, CXCR3, SYK) in normal tissues and tumour tissues from CMM at the protein level. The results showed that compared with normal tissues, these genes were significantly overexpressed in tumour tissues (Fig. 11).

Cutaneous malignant melanoma (CMM) is a skin-related devastating cancer characterized by complicated internetwork, and many genetic factors participate in its progression[1, 5]. For years, the traditional therapeutic methods, such as immunotherapies, chemotherapies, and targeted-therapies, have improved patients’ survival to some extent[17]. However, early-stage diagnostic markers and effective therapies for CMM patients are unavailable at present. Indeed, many patients in clinical suffered drug resistance and allergic by adapting these therapies and result in poor prognosis[3]. Therefore, more patient-friendly therapeutic and valuable prognostic molecules must be established. In recent years, owing to bioinformatics methods evolving that ensure researchers have access to the public database to explore and offer a new angle to carcinoma-involved molecules.

In this study, we divided all the cases downloaded from the TCGA platform into two different groups marked high and low immune/stromal scores, respectively. Then, 908 DEGs (900 upregulation and eight down-regulation) were screened from the immune score and stromal score based on the ESTIMATE algorithm[8]. The functional enrichment analysis indicated that immune response and inflammatory response were firmly involved in CMM tumorigenesis. GO annotation mainly enriched in T cell activation (GO:0042110), regulation of lymphocyte activation (GO:0051249), and regulation of T cell activation (GO:0050863). We conclude that immune cells in extracellular space are of significant importance to tumorigenesis. The results of KEGG pathways in DEGs, including Cytokine-cytokine receptor interaction (hsa04060), Hematopoietic cell lineage (hsa04640) and Cell adhesion molecules (CAMs) (hsa04514), which indicated cytokine and its receptor, coupled with immune cells, played a dominant role in the progression of CMM. Then, the PPI network unveiled the interaction network among these DEGs and demonstrated the intricate interaction of them. We identified ten hub genes and three major modules based on scores, respectively. The core genes are as follows: CCL5, CXCL10, CD74, CXCL9, CD247, IL10, CXCL11, CXCR3, CD3E, and SYK.

C–X–C motif ligand 9,10,11 (CXCL9, CXCL10, CXCL11) is a member of the CXC family chemokine. They are principally expressed by immune cells such as macrophages and T cells, and they play a significant role in immune cell infiltration into the tumour bed. They are supposed to act as a tumour-suppressive molecule[18, 19]. CXCR3, the receptor for CXCL9 and CXCL10, is highly expressed on CD8⁺ and CD4⁺ immune-infiltrating lymphocytes, it is highly likely that local production of CXCL9/10 can modulate T-cell recruitment and activation in human cancers[18–20]. These chemokines are increased following expose to anti-PD-1 and that they are crucial for therapeutic activity and correlate with patient prognosis[21]. The expression of CXCL11 correlated with improved responses to atezolizumab and was the most significantly upregulated gene in macrophages responding to immune checkpoint blockade and its function highly relative to CXCL9 and CXCL10[22, 23]. C–C motif ligand 5, CCL5, is a small protein that belongs to a large family of cytokines and displays chemotaxis activity for it is involved in promoting the migration of several leukocytes into inflammation sites. CCL5 is secreted by a wide variety of cells, including T cells, NKs, and some tumour cells. Previous data showed that targeting CCL5 was sufficient to inhibit the infiltration of NK cells significantly and subsequently enhance the tumour growth[24, 25]. Macrophage migration inhibitory factor (MIF) is an inflammatory cytokine and serves as a regulator of the innate immune system. Studies have shown that MIF induces an immune-suppressive environment that supports melanoma progression and metastasis. CD74 is the central receptor for MIF and the invariant chain of the MHC class II, which plays an essential role in antigen presentation. MIF and CD74 are attractive targets for immunotherapy. MIF binding to CD74 activates the PI3K/AKT and MAPK signalling pathways, and both these pathways have been related to monocyte immune-suppression and macrophage M2-like polarization to regulate the immune response against metastatic melanoma[26–28]. Spleen tyrosine kinase (SYK) functions as a non-receptor kinase, mediating signal transduction of cellular transmembrane receptors and act as immunoreceptors and integrins. SYK has been demonstrated to be a critical regulator of the target of rapamycin (mTOR) activity in B-cell lymphoma[29, 30]. IL-10 is a potent anti-inflammatory molecule produced by innate and adaptive immune cells, including T cells, NK cells, as well as tumour cells. IL-10 is an immune-regulatory cytokine that may exert immune-stimulatory effects on CD8⁺ T cells, depending on their state of activation[21]. CD247 plays a crucial role in triggering the TCR signalling pathway, and high expression levels of CD247 was associated with poor overall survival in lower-grade glioma[2]. Nevertheless, the correlation of CD247 and CMM patients’ overall survival has been no reported in recent years.

Then, we used a univariate Cox model to determine prognostic-related genes highly correlated with overall survival (OS). We combined the identified prognostic-related genes with the immune-related DEGs from the PPI network, and CD247 was the intersection gene for further explored, pair-difference analysis between tumour cases and normal cases and overall survival curve was performed for CD247. The results showed that CD247 highly expressed in tumour samples and possessed a better prognosis than lower expressed samples. In these findings, WGCNA analysis showed that 11 modules were identified. Then, the relationship between each module traits was studied. Specifically, genes belong to module brown had the highest positive correlation with CMM patients’ prognosis, and CD247 included in this module, which suggested CD247 higher expressed samples have a better prognosis.

Moreover, CMM patients with high expression of the ten genes (CCL5, CXCL10, CD74, CXCL9, CD247, IL10, CXCL11, CXCR3, CD3E, and SYK) were associated with better OS. IHC results from the HPA verified that IL10, CXCR3, SYK, CD247, CD74, CXCL11 were significantly overexpressed in CMM tissues. The results suggested that CD247 was an essential biomarker in CMM tumorigenesis, progress, and prognosis. Also, the expressions of CD247 was positively correlated with six immune cells infiltration (CD8⁺ T cell, MacrophagesM1, activated CD4 memory T cell, regulatory T cell, helper follicular T cell, and activated NK cell). Previous studies suggested that high levels of immune cell infiltration are associated with favourable outcomes. Increasing levels of CD8⁺ T cells, MacrophagesM1, activated CD4 memory T cell, NK cell, were related to longer survival time in CMM patients. More critical, CMap analysis tool was used to search for potential drugs or compounds to cure CMM. It was reported that L-securinine inhibits prostate cancer cell growth and metastasis via regulating mitochondrial and AGTR1/MEK/ERK apoptotic pathways and may be a promising candidate medicine for CMM[31]. However, there is no relevant research on the effect of these compounds on CMM.

Even though the paramount role of TME in CMM has been reported, many engaged in the immune infiltrating CD8⁺ T cell. Notably, our study embraces several merits. First, we examined the correlation between immune infiltrating macrophages (M0, M1, M2) and CD247 expression, which is indicated that macrophages, an indispensable part of the microenvironment (TME), was involved in CMM tumorigenesis and progression. Second, the RNA-seq profiling matrix was collected from the TCGA platform, which contained 468 CMM cases. We explored the OncoLnc website, the HPA resource, and the GEPIA database to confirm the results. Finally, the TME and immune-related molecules interact in a changing procession, which subjected to numerous factors, such as immune cells and genetic aspects. We thoroughly investigated the underlying association of TME and immune-linked genes from the element of CMM microenvironment, stromal/immune condition, CD247, and other prognostic-relevant genes and immune cells infiltrating.

In summary, our study mainly illustrated the interaction of TME and immune-linked genes. We finally labelled CD247 as a promising prognostic biomarker for cutaneous cancer.

CMM: Cutaneous malignant melanoma; DEGs: Differentially expressed genes; MEs: Module eigengenes; WGCNA: Weighted gene co-expression network analysis; MCODE: Molecular complex detection; TME: Tumor microenvironment; TCGA: The Cancer Genome Atlas; KEGG: The Kyoto Encyclopedia of Genes and Genomes pathways; PPI: Protein-Protein Interaction; OS: overall survival; GEPIA: the Gene Expression Profiling Interactive Analysis; HPA: the Human Protein Atlas; ECM: extracellular matrix molecules; TOM: topological overlap matrix; STRING: the Search Tool for the Retrieval of Interacting Genes; CMap: ConnectiveMap;

Ethics approval and consent to participate

Not applicable. All data in this study are publicly available.

Consent for publication

Not applicable.

Availability of data and materials

All analyzed data are included in this article and its supplementary materials. The original data are available upon request to the corresponding author.

Competing interests

The authors declare that they have no competing interests.

Funding

The present study was supported by the grants from National Natural Science Foundation of China (NSFC 81570044), Shandong Provincial Natural Science Foundation (2015ZRE27075), and Shenzhen Basic Research Program (No. JCYJ20190807093801661) for Huan-Liang Wang. The funds mentioned were used for

the design of the study, the collection, analysis, and interpretation of the data as well as the writing of the manuscript.

Authors’ contributions

YM designed the research. SSC, YYF and RYM collected and processed the data. YM and HLW wrote the manuscript. All the authors have read and approved the final manuscript.

Acknowledgements

We thank Dr. Sheng-Qiang Li for his kind suggests.

Author details

¹Department of Anesthesiology, Qilu Hospital of Shandong University, Jinan, 250012, China.

²Shenzhen Research Institute of Shandong University, Shenzhen, 518058, China.

Fang S, Xu T, Xiong M, Zhou X, Wang Y, Haydu LE, Ross MI, Gershenwald JE, Prieto VG, Cormier JN et al: Role of Immune Response, Inflammation, and Tumor Immune Response-Related Cytokines/Chemokines in Melanoma Progression. J Invest Dermatol 2019, 139(11):2352-2358 e2353.
Wang Q, Li P, Wu W: A systematic analysis of immune genes and overall survival in cancer patients. BMC cancer 2019, 19(1):1225.
Weiss SA, Han SW, Lui K, Tchack J, Shapiro R, Berman R, Zhong J, Krogsgaard M, Osman I, Darvishian F: Immunologic heterogeneity of tumor-infiltrating lymphocyte composition in primary melanoma. Hum Pathol 2016, 57:116-125.
Ishii G, Ochiai A, Neri S: Phenotypic and functional heterogeneity of cancer-associated fibroblast within the tumor microenvironment. Adv Drug Deliv Rev 2016, 99(Pt B):186-196.
Casey SC, Amedei A, Aquilano K, Azmi AS, Benencia F, Bhakta D, Bilsland AE, Boosani CS, Chen S, Ciriolo MR et al: Cancer prevention and therapy through the modulation of the tumor microenvironment. Semin Cancer Biol 2015, 35 Suppl:S199-S223.
Wang C, Wang X, Liu J, Huang Z, Li C, Liu Y, Sang X, Yang L, Wang S, Su Y et al: Embryonic stem cell microenvironment suppresses the malignancy of cutaneous melanoma cells by down-regulating PI3K/AKT pathway. Cancer Med 2019, 8(9):4265-4277.
Colaprico A, Silva TC, Olsen C, Garofano L, Cava C, Garolini D, Sabedot TS, Malta TM, Pagnotta SM, Castiglioni I et al: TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA data. Nucleic Acids Res 2016, 44(8):e71.
Ke ZB, Wu YP, Huang P, Hou J, Chen YH, Dong RN, Lin F, Wei Y, Xue XY, Ng CF et al: Identification of novel genes in testicular cancer microenvironment based on ESTIMATE algorithm-derived immune scores. J Cell Physiol 2020.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK: limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015, 43(7):e47.
Niemira M, Collin F, Szalkowska A, Bielska A, Chwialkowska K, Reszec J, Niklinski J, Kwasniewski M, Kretowski A: Molecular Signature of Subtypes of Non-Small-Cell Lung Cancer by Large-Scale Transcriptional Profiling: Identification of Key Modules and Genes by Weighted Gene Co-Expression Network Analysis (WGCNA). Cancers (Basel) 2019, 12(1).
Araujo FA, Barh D, Silva A, Guimaraes L, Ramos RTJ: GO FEAT: a rapid web-based functional annotation tool for genomic and transcriptomic data. Sci Rep 2018, 8(1):1794.
Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K: KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res 2017, 45(D1):D353-D361.
Otasek D, Morris JH, Boucas J, Pico AR, Demchak B: Cytoscape Automation: empowering workflow-based network analysis. Genome Biol 2019, 20(1):185.
Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, Nair VS, Xu Y, Khuong A, Hoang CD et al: The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med 2015, 21(8):938-945.
Deng JL, Xu YH, Wang G: Identification of Potential Crucial Genes and Key Pathways in Breast Cancer Using Bioinformatic Analysis. Front Genet 2019, 10:695.
Vella D, Marini S, Vitali F, Di Silvestre D, Mauri G, Bellazzi R: MTGO: PPI Network Analysis Via Topological and Functional Module Identification. Sci Rep 2018, 8(1):5499.
Sadozai H, Gruber T, Hunger RE, Schenk M: Recent Successes and Future Directions in Immunotherapy of Cutaneous Melanoma. Front Immunol 2017, 8:1617.
Doron H, Amer M, Ershaid N, Blazquez R, Shani O, Lahav TG, Cohen N, Adler O, Hakim Z, Pozzi S et al: Inflammatory Activation of Astrocytes Facilitates Melanoma Brain Tropism via the CXCL10-CXCR3 Signaling Axis. Cell Rep 2019, 28(7):1785-1798 e1786.
House IG, Savas P, Lai J, Chen AXY, Oliver AJ, Teo ZL, Todd KL, Henderson MA, Giuffrida L, Petley EV et al: Macrophage-Derived CXCL9 and CXCL10 Are Required for Antitumor Immune Responses Following Immune Checkpoint Blockade. Clin Cancer Res 2020, 26(2):487-504.
Harlin H, Meng Y, Peterson AC, Zha Y, Tretiakova M, Slingluff C, McKee M, Gajewski TF: Chemokine expression in melanoma metastases associated with CD8+ T-cell recruitment. Cancer Res 2009, 69(7):3077-3085.
Sun Z, Fourcade J, Pagliano O, Chauvin JM, Sander C, Kirkwood JM, Zarour HM: IL10 and PD-1 Cooperate to Limit the Activity of Tumor-Specific CD8+ T Cells. Cancer Res 2015, 75(8):1635-1644.
Hong M, Puaux AL, Huang C, Loumagne L, Tow C, Mackay C, Kato M, Prevost-Blondel A, Avril MF, Nardin A et al: Chemotherapy induces intratumoral expression of chemokines in cutaneous melanoma, favoring T-cell infiltration and tumor control. Cancer Res 2011, 71(22):6997-7009.
Mauldin IS, Wages NA, Stowman AM, Wang E, Smolkin ME, Olson WC, Deacon DH, Smith KT, Galeassi NV, Chianese-Bullock KA et al: Intratumoral interferon-gamma increases chemokine production but fails to increase T cell infiltration of human melanoma metastases. Cancer Immunol Immunother 2016, 65(10):1189-1199.
Mgrditchian T, Arakelian T, Paggetti J, Noman MZ, Viry E, Moussay E, Van Moer K, Kreis S, Guerin C, Buart S et al: Targeting autophagy inhibits melanoma growth by enhancing NK cells infiltration in a CCL5-dependent manner. Proc Natl Acad Sci U S A 2017, 114(44):E9271-E9279.
Voshtani R, Song M, Wang H, Li X, Zhang W, Tavallaie MS, Yan W, Sun J, Wei F, Ma X: Progranulin promotes melanoma progression by inhibiting natural killer cell recruitment to the tumor microenvironment. Cancer Lett 2019, 465:24-35.
Figueiredo CR, Azevedo RA, Mousdell S, Resende-Lara PT, Ireland L, Santos A, Girola N, Cunha R, Schmid MC, Polonelli L et al: Blockade of MIF-CD74 Signalling on Macrophages and Dendritic Cells Restores the Antitumour Immune Response Against Metastatic Melanoma. Front Immunol 2018, 9:1132.
Halaby MJ, Hezaveh K, Lamorte S, Ciudad MT, Kloetgen A, MacLeod BL, Guo M, Chakravarthy A, Medina TDS, Ugel S et al: GCN2 drives macrophage and MDSC function and immunosuppression in the tumor microenvironment. Sci Immunol 2019, 4(42).
Imaoka M, Tanese K, Masugi Y, Hayashi M, Sakamoto M: Macrophage migration inhibitory factor-CD74 interaction regulates the expression of programmed cell death ligand 1 in melanoma cells. Cancer Sci 2019, 110(7):2273-2283.
Griss J, Bauer W, Wagner C, Simon M, Chen M, Grabmeier-Pfistershammer K, Maurer-Granofszky M, Roka F, Penz T, Bock C et al: B cells sustain inflammation and predict response to immune checkpoint blockade in human melanoma. Nat Commun 2019, 10(1):4186.
Tang H, Xu X, Xiao W, Liao Y, Xiao X, Li L, Li K, Jia X, Feng H: Silencing of microRNA-27a facilitates autophagy and apoptosis of melanoma cells through the activation of the SYK-dependent mTOR signaling pathway. J Cell Biochem 2019, 120(8):13262-13274.
Zhang D, Liu H, Yang B, Hu J, Cheng Y: L-securinine inhibits cell growth and metastasis of human androgen-independent prostate cancer DU145 cells via regulating mitochondrial and AGTR1/MEK/ERK/STAT3/PAX2 apoptotic pathways. Biosci Rep 2019, 39(5).

Download PDF

Version 1

posted

You are reading this latest preprint version

CD247 Functions as a Prognostic Biomarker for Cutaneous Malignant Melanoma Based on the Analysis of Tumor-immune Microenvironment

Status:

Version 1

Abstract

Figures

Background

Methods

2.1 | Raw data collection

2.2 | Differential analysis of expressed genes

2.3 | WGCNA construction

2.4 | Functional Enrichment annotation and PPI network

2.5 | Clinical correlation analysis of CD247

2.6 | Immune cells infiltration analysis

2.7 | ConnectiveMap (CMap) and overall survival analysis

2.8 | Statistical methods

Results

3.1 | Stromal and immune scores are highly linked with the TNM staging system and prognosis

3.2 | Analysis of Differential expression genes (DEGs) from CMM cases based on stromal and immune scores

3.3 | Weighted gene co-expression network analysis

3.4 | Functional Enrichment annotation and PPI network analysis of DEGs

3.5 | CD247 was determined as the most valuable prognostic gene.

3.6 | The relationship between immune cells infiltration abundance and prognosis of CMM

3.7 | The survival curve of prognostic-related hub genes and ConnectiveMap (CMap) analysis

3.8 | Mining the OncoLnc database and HPA platform to validate obtained prognostic genes

Discussion

Conclusion

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Authors’ contributions

Acknowledgements

Author details

References

Supplementary Files

Status:

Version 1