Development and Validation of a Combined Glycolysis and Immune Prognostic Model for Melanoma

Background Glycolytic effects and immune microenvironments play important roles in the development of melanoma. However, reliable biomarkers for prognostic prediction of melanoma as based on glycolysis and immune status remain to be identified. Methods Glycolysis-related genes (GRGs) were obtained from the Molecular Signatures database and immune-related genes (IRGs) were downloaded from the ImmPort dataset. Prognostic GRGs and IRGs in the TCGA (The Cancer Genome Atlas) and GSE65904 datasets were identified. Least absolute shrinkage and selection operator (LASSO) Cox regression and multivariate Cox regression were used for model construction. Glycolysis expression profiles and the infiltration of immune cells were analyzed and compared. Finally, in vitro experiments were performed to assess the expression and function of these CIGI genes. Results Four prognostic glycolysis- and immune-related signatures (SEMA4D, IFITM1, KIF20A and GPR87) were identified for use in constructing a comprehensive glycolysis and immune (CIGI) model. CIGI proved to be a stable, predictive method as determined from different datasets and subgroups of patients and served as an independent prognostic factor for melanoma patients. In addition, patients in the high-CIGI group showed increased levels of glycolytic gene expressions and exhibited immune-suppressive features. Finally, SEMA4D and IFITM1 may function as tumor suppressor genes, while KIF20A and GPR87 may function as oncogenes in melanoma as revealed from results of in vitro experiments. Conclusion In this report we present our findings on the development and validation of a novel prognostic classifier for use in patients with melanoma as based on glycolysis and immune expression profiles.


INTRODUCTION
Skin cutaneous melanoma (SKCM) is a cancer resulting from malignant transformations of melanocytes in the basal layer of the skin. This condition represents one of the most aggressive types of skin cancer worldwide, with factors such as ultraviolet light (UV) and malignant transformations of moles being high-risk factors for SKCM (1). At present, most treatments for SKCM involve surgical resection removal, radiotherapy and chemotherapy (2). Although significant progress has been achieved with these treatments, more than 95% of metastatic melanoma patients die within one year (3). Therefore, there is an urgent need to identify prognostic biomarkers, which would provide clinicians the ability to promptly and accurately predict the clinical outcome of melanoma as well as initiate a protocol for a personalized treatment regimen.
Abnormal immune microenvironments and tumor metabolic reprogramming represent two significant characteristics of tumors (4). Tumor metabolic reprogramming is one of the critical mechanisms involved with regulating the immune microenvironment (5,6). Tumor metabolic reprogramming includes enhanced aerobic glycolysis (Warburg effect), increased glucose uptake and consumption, enhanced lipid and protein synthesis and enhanced glutamine uptake and catabolism, with glycolysis as the main feature (7)(8)(9). Even with sufficient oxygen availability, tumor cells continue to metabolize glucose mainly through glycolysis to produce lactic acid (10,11). The capacity for this glycolysis within tumor cells along with its production of lactic acid to participate in the regulation of immune microenvironments have become notable areas of research interest in recent years. It has been reported that lactic acid, as secreted from tumor cells, can inhibit the cytolytic ability of CD8+ effector T cells, but it may also be used by Treg cells to support cell metabolism (12). The lactic acid produced by melanoma cells can reduce the immune surveillance ability of T and NK cells by inhibiting the nuclear factor produced by the activated T cell (NFAT) dependent IFN-g (13). In addition, the findings that glycolysis is the main metabolic pathway required after T cell activation, and enhances the eradicating effects of T cells (14,15), indicates that glycolytic effects play different roles in immune versus tumor cells. Although a clear connection exists between glycolytic effects and immune microenvironments, scant attention has been directed toward examining this relationship in any detail.
In this study, we constructed and validated a comprehensive index of glycolysis and immune (CIGI) model as based on glycolysis-(GRGs) and immune-(IRGs) related genes. This CIGI model shows stable, prognostic prediction performance in different data sets and in different subgroups of melanoma patients. In additional, we demonstrate that CIGI was correlated with glycolysis and immune status. Finally, the expression and function of CIGI genes in melanoma were evaluated in in vitro experiments.

Data Collection
Gene expression profile and clinical follow-up information of RNA sequencing samples from patients with SKCM were downloaded from the TCGA database, and the GSE65904 cohort was downloaded from the Gene Expression Omnibus (GEO) database. The RNA sequencing data from TCGA-SKCM were preprocessed as follows: 1) the samples without follow-up information were removed; 2) ENSEMBL ID were converted to gene symbols; 3) Remove the samples with OS <30 days; 4) Genes with expression level is lower than 1 and the proportion is higher than 50% in all samples were eliminated. The GEO cohort was processed via the following steps: 1) the samples without clinical follow-up information were removed; 2) ENSEMBL IDs were converted to gene symbols; 3) the probes corresponding to multiple genes were removed; 4) the median of multiple gene expression values was used. TCGA database and GSE65904 dataset were selected as they contained the largest sample set in the same platform with detailed followup information on melanoma. Glycolysis-related genes were identified in the Molecular Signatures database, while immunerelated genes were downloaded from the ImmPort dataset.

Construction and Validation of the CIGI
The independent prognostic predictors among glycolysis-and immune-related genes were selected via univariate analysis based on the "survival" package. LASSO analysis and the stepwise Cox proportional hazards regression model were used to construct CIGI. High-and low-CIGI groups were differentiated as based on the optimized risk value. Kaplan-Meier survival analyses were used to analyze differences in overall survival between the highand low-CIGI groups. Time-dependent ROC curve analysis was used to evaluate the predictive value of the CIGI. Univariate and multivariate Cox regression analyses were performed to explore the independent prognostic value of the CIGI.

Potential Regulatory Pathways Analysis
Single-sample gene set enrichment analysis (ssGSEA) based on the "GSVA" package was used to quantify the scores of pathways in each sample. And, the ssGSEA score from each sample was used to analyze the relationship between CIGI and potentially regulatory pathways.

Immune Infiltration Analysis
The infiltration levels of 28 immune cell types in the TCGA dataset were assessed using ssGSEA analysis. The gene markers of 28 immune cells were obtained from a previous study (16).

Univariate and Multivariate Cox Analysis
Univariate and multivariate COX regressions were used to analyze the relationship between CIGI and other variables and the clinical prognosis of patients. (1:2000; ab233545; Abcam), KIF20A (1:100; ab70791; Abcam) and GPR87 (1:100; ab188901; Abcam). Each sample was independently assessed by two pathologists and scored using a semiquantitative scoring system with histoscores ranging from 0 (minimum) to 300 (maximum).

IHC (Immunohistochemistry) Analysis
This, as well as any other portions of this study involving human tissue samples, were approved by the Ethics Committee of the First Hospital of China Medical University.

Cell Culture
The human melanocytes PIG1 cell line, melanoma cell line, A375, A875, and MeWo were obtained from the China Infrastructure of Cell Line Resource and was cultured in DMEM (10% FBS) at 37°C in a humidified 5% CO 2 incubator.

Colony-Forming Experiments
In order to examine the effects of SEMA4D, IFITM1, KIF20A and GPR87 expression on human melanoma cell proliferation, A375 cells (500/well) transfected with NC-siRNA or siRNAs were added to the 12-well plates. After two weeks, the number of colonies were counted. siRNA, GPR87-specific siRNA). The culture medium and cells were collected after 48 h. Lactate levels in the medium were determined with use of the lactate assay kit (ab65331, Abcam), ATP levels with the ATP assay kit (ab83355, Abcam) and glucose uptake levels with the glucose uptake assay kit (ab136955, Abcam). All determinations were normalized with cell numbers.

Extracellular Acidification Rate (ECAR)
Measurement of ECAR was performed using the Seahorse XF96 Flux Analyser (Seahorse Bioscience) following instructions provided by the manufacturer. In short, after transfection with NC-siRNA or siRNAs (KIF20A-specific siRNA, GPR87-specific siRNA), A375 cells (1.5 × 10 4 /well) were cultured in 96-well plates. The culture medium was replaced with the test buffer prior to the ECAR test. Cells were then incubated with buffered medium containing 10mM glucose, 1mM oligomycin and 50mM 2-DG.

Statistical Analysis
The Kaplan-Meier analysis was used to compare differences in OS between the high-and low-CIGI groups. Statistical comparisons between two groups were performed using the Student's two-tailed t-test. A p-value of < 0.05 was required for results to be considered as statistically significant.

Development of CIGI in Melanoma
TCGA and GSE65904 datasets were included in this study. Univariate Cox regression analysis was conducted to select genes significantly related with prognosis based on these TCGA and GSE65904 datasets. A total of 6155 genes were significantly related to the overall survival of melanoma patients in the TCGA dataset, with hazard ratios of the top 10 genes(sorted by p value) presented in Figure 1A. A total of 2835 genes were significantly related to the overall survival of melanoma patients in the GSE65904 dataset, with hazard ratios of the top 10 genes(sorted by p value) presented in Figure 1B. A Venn diagram indicating that 182 prognostic GRGs and IRGs were identified in the GSE14520 and TCGA cohorts as shown in Figure 1C. As an excessively large number of genes is not conducive to clinical detection, we further narrowed the range of genes. We employed the Lasso regression analysis and the resulting change trajectory of each independent variable is shown in Figure 1D. With the gradual increase of lambda, the number of independent variable coefficients gradually increased to zero. Five-fold cross-validation was used to build the model, and the confidence interval under each lambda. Thus, we selected 17 genes at lambda= 0.0748 as the candidate genes. 17 genes were as follow: SEMA4D, PDGFRB, C5, PSMC6, CCL8, CNTFR, IL27RA, PIK3R2, IFITM1, KIF20A, GPR87, LEP, MAP2K1, COL5A1, SSTR2, SEMA6A and KLRD1. To optimize this model and identify only the most predictive genes, a stepwise Cox proportional hazards regression model was used, which resulted in a final set of 4 genes ( Figure 1E). As a result of these final analyses the CIGI was constructed: CIGI = (-0.164 × SEMA4D expression) + (-0.15 × IFITM1 expression) + (0.278 × KIF20A expression) + (0.137 × GPR87 expression).

Prognostic Analysis of CIGI in the TCGA Dataset
The risk score was calculated and distribution of the TCGA cohort is shown in Figure 2A. As based on the optimized risk value, patients were assigned to either a high-or low-risk CIGI group. Results from the Kaplan-Meier survival analyses indicated that the overall survival of patients in the high-CIGI group was significantly lower than that of patients in the low-risk group ( Figure 2B; p < 0.0001). The ROC curve analyses showed that the AUC values for the 1-, 3-, and 5-year survival rates were 0.738, 0.655 and 0.702, respectively ( Figure 2C). Additionally, patients with low-risk scores experienced a significantly longer diseasefree survival (DSS, p <0.001; Figure 2D) and progression-free interval (PFI, p = 0.0015; Figure 2E).

Verification of CIGI Based on the GSE65904 Dataset
Next, the stability and reliability of CIGI in the GSE65904 dataset were assessed. The CIGI distribution in the GSE65904 cohort is shown in Figure 3A. Results from further prognostic analysis using Kaplan-Meier revealed that the overall survival of melanoma patients in the high-CIGI group was significantly lower than that in the low-CIGI group (p=0.00054; Figure 3B). ROC curve analyses showed that AUC values for the 1-, 3-and 5-year survival rates were 0.622, 0.631 and 0.58, respectively ( Figure 3C).

Prognostic Value of CIGI in Different Subgroups of Patients With Melanoma
To further evaluate the clinical application value of CIGI, the prognostic value of CIGI in melanoma patients with different clinical characteristics was analyzed. As summarized in Figure 4, among melanoma patients with different clinical characteristics the OS of the high-CIGI group was significantly lower than that of the low-CIGI group (Figures 4A-L) Figure 5B). Taken together, these results suggest that a high CIGI value served as an independent factor for prognosis in melanoma patients.

Nomogram Construction
A nomogram was generated to provide a quantitative method for predicting the probability of 1-, 3-, and 5-year OS in patients with melanoma, which could then be used in clinical practice. The nomogram integrated clinicopathological features and CIGI as based on results of the multivariate Cox regression analysis ( Figure 6A). The c-index value of the nomogram was 0.74, indicating a satisfactory overlap with actual observations. Calibration curves showing calibration points in 5-years demonstrate a high degree of coincidence with the standard  curves, indicating that the model provides an effective level of prediction performance ( Figure 6B). In additional, results from the decision curve analysis indicate that the nomogram can better predict OS than that obtained using single clinicopathological features ( Figure 6C). Therefore, this nomogram, as based on CIGI, could be used to predict the prognosis of melanoma patients in clinical practice.

Glycolysis Profile in the CIGI
Next, the correlation between CIGI and glycolysis was assessed. As reported previously (17), multiple genes are involved in encoding glucose transporters and critical kinases (glycolytic genes) that directly regulate cell glycolysis ( Figure 7A). Here, we analyzed the changes in mRNA of these genes between the lowversus high-CIGI group in the TCGA and GSE65904 databases. With the exception of HK2, PFKL and SLC2A1, the high-CIGI group showed higher expressions of other glycolytic genes than that of the low-CIGI group as based on the TCGA dataset ( Figure 7B). Similarly, within the GSE65904 dataset, the expression of glycolytic genes were again higher in the high-CIGI group, with the exception of PFKL, PGAM1 and SLC2A1 ( Figure 7C). In additional, ssGSEA scores of the "HALLMARK GLYCOLYSIS" pathway via GSVA analysis, as based on the TCGA and GSE65904 datasets were also generated ( Figures 7D, E). From this analysis, we found that activity of the glycolysis pathway increased as a function of increases in CIGI. Overall, these results suggest a glycolysis-overexpression status is present in the high-CIGI group.

Immune Profile in the CIGI
In this experiment, the relationship between CIGI and the immune status of melanoma was evaluated. To accomplish this goal a correlation analysis was first performed between CIGI and immune infiltrate cells, with the result that infiltration levels of 28 immune cell types were obtained using ssGSEA. With the exclusion of CD56bright natural killer and Type 2 T helper cells, the remaining immune cells were all significantly down-regulated in the high-CIGI group ( Figure 8A). Next, the relationship between CIGI and tumor immune microenvironments were assessed. Stromal and immune scores were used to estimate tumor immune microenvironments. As shown in Figures 8B, C, CIGI was negatively correlated with immune (R=−0.46, p<0.001) and stromal (R=−0.36, p<0.001) scores. Finally, the relationship between CIGI and immune-related pathways were determined. With this analysis, we found that the activity of the Interferon-a response, Interferon-g response, TNF-a via NF-KB, Inflammatory response and the IL2-STAT5 signaling pathway were all increased with increases in CIGI ( Figure 8D). Overall, the results of these analyses indicate that the high-CIGI group exhibited immunesuppressive features.

Protein Level Validation and Functional Analysis of CIGI Genes
In Figures 9A, B the mRNA and protein level of SEMA4D and IFITM1 were decreased in A375, A875, and MeWo cells compared with melanocytes PIG1 cells. And, the mRNA and protein level of KIF20A and GPR87 were upregulated in A375, A875, and MeWo cells compared with melanocytes PIG1 cells.
To verify the protein expressions of SEMA4D, IFITM1, KIF20A and GPR87 in melanoma tissue, 20 melanoma and paired normal tissue samples were compared. Results from the immunohistochemistry assay revealed that SEMA4D and IFITM1 were down-regulated, while KIF20A and GPR87 were over-expressed in melanoma tissue ( Figure 9C). When analyzing the potential biological functions of these proteins in melanoma, we found that silencing SEMA4D or IFITM1 promoted, while silencing KIF20A or GPR87 inhibited the proliferation of melanoma cells in vitro (Figures 9D, E). These results suggest that SEMA4D and IFITM1 may function as tumor suppressor genes while KIF20A and GPR87 may function as oncogenes in melanoma.
KIF20A and GPR87 Regulate the Glycolysis Ability of Melanoma Cells As KIF20A and GPR87 are glycolysis-related genes with poor prognostic potential for patients with melanoma, we next examined the effects of KIF20A and GPR87 on glycolytic ability within melanoma cells. After transfection with KIF20A-siRNA or GPR87-siRNA, levels of lactate, ATP and glucose uptake as well as extracellular acidification rates were all significantly reduced in these melanoma cells (Figures 10A-D). Based on these results, it seems likely KIF20A and GPR87 may function as oncogenes through their capacity to regulate glycolysis within melanoma cells. A large number of studies have shown that the AKT/LDHA pathway plays a vital role in the process of cellular glycolysis (18)(19)(20). We tried to explore the regulatory effects of KIF20A and GPR87 on the AKT/LDHA pathway. The results showed that after knocking out KIF20A and GPR87, the AKT/LDHA pathway activity of A375 cells was significantly inhibited ( Figure 10E). These results suggest that: KIF20A and GPR87 may regulate the glycolytic ability of melanoma cells by regulating the AKT/LDHA pathway.

DISCUSSION
Melanoma is a highly aggressive cancer, with a poor prognosis largely due to limited treatments currently available (21). Therefore, prognosis signatures for melanoma patients are sorely needed. With the development of bioinformatics and sequencing technology, several potential prognostic evaluation programs for melanoma patients have been generated (22)(23)(24). However, most of the parameters analyzed in these programs originate from the genome or transcriptome, with no consideration of biological processes. Therefore, these models do not adequately reflect the disease characteristics of melanoma. Glycolysis and immune microenvironments represent two significant biological hallmarks of tumors that have been demonstrated to be of value in evaluating the prognosis of patients with melanoma (25,26). In this study, we included GRGs and IRGs to construct a CIGI as based on gene expression data in public databases. As demonstrated within different data sets and subgroups of patients, CIGI exhibits an effective degree of predictive performance, and can serve as an independent prognostic factor for melanoma patients. In this way, the accuracy displayed by CIGI for the prognosis of melanoma indicates a great potential for clinical application. Several studies have constructed prognostic models of melanoma based on gene expression profiles from public databases. Wang et al. identified an 8-genes risk signature in melanoma patients based on entire genome (24). Song et al. constructed a 12-gene signature for survival prediction in malignant melanoma patients (27). However, these studies were derived from the entire genome or transcriptome with no consideration of biological processes. Consequently, they were simply mathematic models that did not reflect the intrinsic character of the cancer itself. Recently, researchers try to build a prognostic risk model based on tumor cell characteristics such as immune, metabolism, hypoxia, etc. Xue et al. identified an immunerelated signature for predicting prognosis in melanoma, which contained 23 immune-related gene pairs (28). Shou et al. determined of hypoxia signature (20 hypoxia-related genes) to predict prognosis in melanoma (29). Our research is the first to construct a prognostic model for the two major tumor characteristics of glycolysis and immunity, which can simultaneously reflect the tumor's immune status and glycolysis changes. In addition, our model contains less genes, which is convenient for clinical application. Importantly, the function of genes in our model has been confirmed based on experiments.
Results from increasing numbers of studies have indicated that metabolic changes in tumor microenvironments can inhibit immune cell infiltrations and other anti-tumor immunity processes through the production of immunosuppressive metabolites (30,31). Glycolysis represents the main feature of tumor metabolism. Glycolysis forms an acidic tumor microenvironment by increasing lactic acid efflux, inhibiting anti-tumor responses mediated by T cells and inhibiting the activity of tumor infiltrating myeloid cells (12,13). Our current findings reveal that melanoma patients in the CIGI-high group show higher expressions of glycolytic genes and exhibit immunesuppressive features. The CIGI we constructed contains four genes, SEMA4D, IFITM1, KIF20A and GPR87. Among these, SEMA4D and IFITM1 are immune-related genes, while KIF20A and GPR87 are glycolysis-related genes. SEMA4D functions as a cell surface receptor for PLXNB1 and PLXNB2 and plays an important role in cell-cell signaling (32). SEMA4D has the capacity of inducing Bcells to aggregate and improves their viability (33),and high expression levels of SEMA4D are significantly correlated with poor prognosis for a number of cancers including colon (34), ovarian epithelial (35) and cervical (36) as well as soft tissue sarcoma (37). With regard to IFITM1, the protein encoded by this gene is an interferon-induced antiviral protein, which mediates the innate immunity of cells to influenza A H1N1, West Nile and Dengue viruses by inhibiting their initial replication (38), and plays a key role in the antiproliferative action of IFN-gamma (39). IFITM1 has been reported to be abnormally expressed in tumor tissues and it is an independent prognostic biomarker for patients with acute myeloid leukemia (40), lung adenocarcinoma (41) and gallbladder cancer (42). KIF20A plays a crucial role in cell mitosis and cell migration. Findings from several studies (43)(44)(45) have indicated that KIF20A is abnormally expressed in tumor tissues and associated with poor prognosis for patients with soft tissue sarcoma, ovarian cancer, and breast cancer. Finally, GPR87, which belongs to the G protein-couple receptor family, has been shown to be overexpressed in cancers such as pancreatic (46), non-small-cell lung (47), bladder (48) and hepatocellular (49). In this study, we present the first evidence that SEMA4D, IFITM1, KIF20A and GPR87 may possess a prognostic value for melanoma. In specific, the results of our in vitro experiments suggest that SEMA4D and IFITM1 may function as tumor suppressor genes while KIF20A and GPR87 may function as oncogenes in melanoma. Moreover, our findings provide the first indication that KIF20A and GPR87 can regulate the glycolytic ability of melanoma cells.
Although the results of our current study suggest that CIGI can serve as an effective prognostic tool for use in melanoma Immune-related pathways with differential expressions between the high-versus low-CIGI group. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001. ns, no significance.
patients, the limitations associated with this study indicate the need for additional analyses prior to a clinical application of this protocol. First, as all samples used in our study were obtained retrospectively, the inclusion/analysis of prospective samples will need to be included for verification of our findings. Second, we focused our analyses only as related to the prognostic value and clinical significance of CIGI. The other potential functions of SEMA4D, IFITM1, KIF20A and GPR87 in CIGI will require further investigation with use of additional in vivo and in vitro experiments

CONCLUSION
In conclusion, here we present the findings of our endeavors to develop and validate a novel prognostic classifier for use in patients with melanoma as based on glycolysis and immune expression profiles.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: TCGA GSE65904.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Ethics Committee of the First Hospital of China Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YY and YL made substantial contributions to the conception and design, acquisition, analysis, and interpretation of data. RQ and LZ was involved in drafting the manuscript and revising it critically for important intellectual content. All authors contributed to the article and approved the submitted version.