Prognostic Role of Tumor Mutation Burden Combined With Immune Infiltrates in Skin Cutaneous Melanoma Based on Multi-Omics Analysis

Tumor mutation burden (TMB) and tumor infiltrating lymphocytes have been well-recognized as molecular determinants of immunotherapeutic responsiveness in many types of cancer. However, the relationship between TMB with immune infiltrates and their prognostic role are reported occasionally in skin cutaneous melanoma (SKCM). We obtained the somatic mutation data and transcriptome profiles of 454 SKCM patients from The Cancer Genome Atlas (TCGA) database, and analyzed the mutation profiles using “maftools” package. Correlation analysis revealed that lower TMB levels conferred poor survival outcomes, associated with lower age and advanced pathological stage. Differential analysis was conducted to the genome expression between two TMB groups using “limma” package, and we identified four hub TMB-related immune genes including CNTFR, CRABP2, GAL, and PAEP. We further analyzed the underlying relationships of the copy number variations (CNVs) of four hub genes with immune infiltrates in SKCM microenvironment through TIMER database. The results indicated that diverse forms of CNVs carried by hub genes could commonly inhibit immune infiltrates. Based on the CIBERSORT method, we compared the proportions of 22 immune cells in two TMB groups and assessed their prognostic value. The data revealed that infiltrations levels of regulatory T (Treg) cell and dendritic activated cells in high-TMB group were lower than that in low-TMB group, while M1 and M2 macrophages showed the opposite trend, especially the levels of neutrophil and macrophage correlated positively with prognosis of SKCM. Finally, we constructed a TMB Prognostic Index (TMBPI) to evaluate the predictive accuracy of the four hub TMB-related immune genes. The ROC curve was drawn to assess the predictive accuracy with AUC = 0.664 and higher TMBPI conferred poor survival outcomes, which warranted further investigation and larger samples to validate.


INTRODUCTION
Malignant melanoma is one of the most aggressive cancers, the incidence of which is rising worldwide (1)(2)(3). On the basis of the anatomical location, melanoma is subdivided into three subtypes: skin cutaneous melanoma (SKCM), acral melanoma, and mucosal melanoma (4,5). Skin cutaneous melanoma (SKCM) is the major subtype of melanoma in Caucasians, which accounts for more than 90% (6). Nevertheless, the proportion of SKCM is approximately 20% in Asian populations (7). In recent years, advances in immunotherapy have significantly improved survival outcomes of SKCM patients. The immune checkpoint blockade (ICB) targeting PD-1/PD-L1 and CTLA-4/B7-1 have been approved for the treatment of advanced SKCM by Food and Drug Administration (FDA) (8)(9)(10)(11). However, the overall efficacy rate of PD-1 blockade in SKCM is approximately 26 to 44% (8,12), thus indicating more than 50% of SKCM patients are not suitable for ICB therapy. Therefore, identification and characterization of potential biomarkers and their application in combination with immunotherapy are urgently required.
In recent years, there are many well-recognized molecular determinants of immunotherapeutic responsiveness, including PD-L1 expression on tumor (13), microsatellite instability (14), tumor mutation burden (TMB) (15), neo-antigen load (16), and tumor infiltrating lymphocytes (TILs) (17). A number of clinical trials have explored the relationship between PD-L1 expression and immunotherapeutic response in several types of cancer, such as non-small cell lung cancer (NSCLC), esophageal cancer, and SKCM (18). The results of KEYNOTE-042 (NCT02220894) study indicated that NSCLC patients with PD-L1 tumor proportion score (TPS) of 1% or greater could benefit from pembrolizumab monotherapy (19). Subgroup analysis showed that NSCLC patients with PD-L1 TPS of 50% or greater have a better survival outcome compared with patients with PD-L1 TPS of 1 to 49%. The retrospective analysis revealed that positive association of PD-L1 expression levels with a better prognosis in SKCM patients treated with ICB therapy (20). TMB was calculated as (total count of variants)/(the whole length of exons), in which the detected variants included base substitutions, insertions, or deletions across bases. Many studies demonstrated that higher TMB in tumors is inclined to form more neo-antigens that make tumors harboring higher immunogenicity, and thus lead to improved clinical response to immunotherapy (21). A series of clinical trials showed that high TMB confers high response rate and sustainable response to ICB therapy in many types of cancer (22,23).
Nevertheless, PD-L1 expression and TMB are still not perfect biomarkers for immunotherapy, as responses are observed in PD-L1 negative or low TMB patients. A pooled analysis of two trials showed that significant benefit (i.e. including complete and durable responses) with nivolumab in all PD-L1 subgroups, including PD-L1 negative NSCLC patients (24). This may be attributed to the immune status of tumor microenvironment. ICB therapy exerts anti-tumor effects depending on the involvement of TILs, thus the abundance of TILs could be regarded as one biomarker for predicting the immunotherapeutic efficacy (25). Given the limitation of single biomarker, a prognostic model containing various biomarkers may guide immunotherapy more precisely.
A recent study investigated the correlation between TMB and immune signatures in different types of cancer (26). However, a systemic exploration of the relationship between TMB with immune infiltrates in SKCM remains lacking, so we performed this study to explore the prognostic role of TMB and its potential association with immune infiltrates in SKCM.

Muti-Omics Data Acquisition and Processing
First, we obtained the somatic mutation data of SKCM samples from The Cancer Genome Atlas (TCGA) database using the GDC tool (http://portal.gdc.cancer.gov/). Since the raw SNP data were not publicly available, we selected the "Masked Somatic Mutation" data and processed it through VarScan software. We prepared the Mutation Annotation Format of somatic variants and implemented the "maftools" R package which provides various functions to perform most commonly used analysis in cancer genomics and to create feature-rich customizable visualizations (27). Then, we downloaded the transcriptome profiles with HTSeq-FPKM format of all available SKCM samples compared with normal tissues. Moreover, the clinical and pathological data of SKCM patients were also obtained via the GDC tool, including age, gender, tumor-node-metastases (TNM) stage, and follow-up with vital status. According to the American Joint Committee on Cancer classification (AJCC) on Melanoma classification, T, N, and M stage refer to tumor thickness, regional lymph nodes metastases, and distant metastases. Since all the data in this research were form public databases, there was no ethical conflict to declare.

Calculation of TMB Scores and Prognostic Analysis
TMB was defined as the total count of coding errors of somatic genes, base substitutions, deletions, or insertions across per million bases. In this study, we calculated the mutation frequency with number of variants/the length of exons (38 million) for each ample through Perl scripts based on the JAVA platform. The estimated TMB data of SKCM samples are shown in Table S1. Then, we classified the SKCM samples into high-TMB and low-TMB groups according the median data. Kaplan-Meier analysis with log-rank test was subsequently performed to compare the survival difference between two groups. In addition, we further evaluated the associations of TMB levels with clinical characteristics via Wilcoxon rank-sum test. According to the tumor classification of malignant degree, we divided the SKCM samples into two groups with the median level.

Differentially Expressed Genes and Functional Pathways Analysis
According to the TMB levels, we divided the transcriptome data of SKCM samples into high-and low-TMB groups via R software. We selected the "limma" package to conduct the differentially expressed gene (DEGs) analysis in two groups with | Fold change (FC) | >1 and False Discovery Rate (FDR) <0.05. The heatmap plot was drawn using "pheatmap" package. Then, "org.Hs.eg.db" package was used to get the Entrez ID for each DEG and we conducted the Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis with "clusterProfiler," "ggplot2," and "enrichplot" packages. Besides, gene set enrichment analysis (GSEA) software was obtained from the home website (http:// software.broadinstitute.org/gsea/index.jsp ) and worked on JAVA platform. We used the TMB level as the phenotype and selected the "c2.cp.kegg.v7.0 symbols.gmt gene sets" as the reference gene set. Furthermore, we obtained a list of immune related genes from the Immunology Database and Analysis Portal (Immport) to select the differentially expressed immune genes between the two groups using "VennDiagram" package.

Survival Analysis
We selected the top 25 immune genes with |FC| >1 and FDR <0.05 to further assess the prognostic value of differential immune genes in patients with high-and low-TMB levels. The survival Kaplan-Meier analysis was conducted via a "for cycle" R script to identify the hub immune genes associated with survival outcomes and P value was shown in plot. A P value <0.05 was regarded as the statistical significance.

TIMER Database and CIBERSORT Algorithm
We further evaluated the copy number variations (CNVs) of hub immune genes with immune infiltrates in SKCM based on the "SCNA" module of TIMER database (http://cistrome.shinyapps. io/timer/) (28). The known CNV types of hub genes were shown at the right bottom. The distributions of each immune cell subset at each CNV status in SKCM were presented by box plots and the difference of infiltration level in each category versus normal was compared using two-sided Wilcoxon rank sum test with calculated P value.
Meanwhile, we obtained the transcriptome profiles of SKCM patients in two TMB levels and conducted the normalization process using "limma" package. Then, we put the preparation data into subsequent analysis to evaluate the immune infractions of each sample through the CIBERSORT algorithm (R script v1.03), providing an estimation of the abundances of member cell types in a mixed cell population, using gene expression data. The CIBERSORT was still based on a known reference set, providing a set of gene expression features of 22 leukocyte subtypes-LM22. The differential distributions of immune cells in two TMB levels were shown using "pheatmap" package. The Wilcoxon rank-sum test was conducted to precisely assess the differential abundances of immune infiltrates between two TMB levels, which were shown with P value using "vioplot" package.

Prognostic Analysis of Infiltrating Immune Cells in SKCM
Based on the TIMER database, we further performed the multivariate Cox of immune infiltration cells, which was fitted by function coxph() from R package "survival." The hazard ratio (HR) with 95% confidence interval (95% CI) was calculated. Furthermore, batch Kaplan-Meier analysis was conducted to indicate the differential survival outcomes between different levels of immune infiltrates. A P value <0.05 of log-rank test was regarded as the statistical significance.

Construction of TMB Prognostic index (TMBPI) for Hub Immune Genes
We performed the multivariate Cox regression analysis to obtain the respective coefficients (bi) of 4 hub immune genes. As described previously (29), the TMBPI formula was defined as: TMBPI = Ʃ(bi × Expi) (i = 4). Then, we performed the Receiver Operating Characteristic (ROC) curve to evaluate the predictive value of 4 immune signature in SKCM. Moreover, Kaplan-Meier analysis was conducted to compare the survival difference in two groups, where we divided the SKCM patients into high-and lowrisk groups with the median prognostic index as the threshold.

Statistical Analysis
The Cox regression model was performed based on the "survival" package. "Limma" package was utilized to conduct the normalization and differential analysis. Wilcoxon rank-sum test was a non-parametric statistical hypothesis test mainly used for comparisons between two groups and Kruskal-Wallis test was suitable for two or more categories. All statistical analysis was implemented based on the R software (Version 3.6.3). A P value <0.05 was thought to be statistically significant.

Patient Characteristics
Four hundred and fifty-four patients with SKCM were included in this study and their clinicopathological characteristics are listed in Table 1

Landscape of Mutation Profiles in SKCM
All mutation information of each gene in each sample were exhibited in waterfall plot, where various colors with annotations at the bottom represented the different mutation types ( Figure  1). On the whole, these mutations were further summarized in different groups, in which missense mutation accounts for the most fraction ( Figure S1A), single nucleotide polymorphism occurred more frequently than deletion or insertion ( Figure  S1B), and C>T was the most common of single nucleotide variants in SKCM ( Figure S1C). Moreover, we calculated the number of altered bases for every sample and showed the mutation categories with various colors in box plots (Figures S1D, E). Horizontal histogram revealed the top 10 mutated signature in SKCM with percentages, including TTN (72%), MUC16 (67%), BRAF (51%), DNAH5 (49%), PCLO (44%), LRP1B (38%), ADGRV1 (35%), RP1 (33%), ANK3 (32%), and DNAH7 (32%) ( Figure S1F). The coincident and exclusive associations across mutated genes were shown in Figure S1G, in which green represented the co-occurrence and red represented the mutually exclusive relationships. Although BRAF mutation may be exclusive to other gene mutations, the differences did not reach a significance. Meanwhile, the mutated incidences of other genes were shown in Genecloud plot ( Figure  S2). Furthermore, we evaluated the above gene variants through cBioPortal for Cancer Genomics ( Figure S3). The results verified the reliability and accuracy of mutation profiles in SKCM.

Correlation of TMB to Survival Outcomes and Clinicopathological Features of SKCM
After dividing TMB into two groups according to the median level, we analyzed the prognostic significance of TMB for overall survival (OS). The results indicated that SKCM patients in low-TMB group had a significantly shorter OS than that in high-TMB group (P = 0.006; Figure 2A). Furthermore, we examined the relationship between TMB and clinicopathological features. Statistical analysis showed that higher TMB level correlated with higher age (P = 0.002, Figure 2B), lower AJCC-T stage (P = 0.044, Figure 2D), lower AJCC-N stage (P = 0.029, Figure  2E), and early pathological stage (P = 0.023, Figure 2C). However, no significant association was found between TMB and AJCC-M stage ( Figure 2F).

Comparison of Gene Expression Profiles Between Two TMB Groups
The heatmap showed that the genome expression levels commonly decreased in high-TMB group than that in low-TMB group ( Figure 3A). Differential analysis revealed a total of 224 DEGs with |log FC|>1 and FDR <0.05 (Table S2). We then conducted the GO enrichment analysis and these DEGs were mainly in skin development, epidermis development, and cell differentiation crosstalk ( Figure 3B, Table S3). KEGG pathway analysis indicated that the enrichment of TMB-related signature correlated with multiple cancer-related crosstalk, such as PI3K-AKT signaling pathway, and Wnt signaling pathway ( Figure 3C, Table S4). In addition, we also selected partial GSEA results of the top TMB-related items, including cell cycle, pyrimidine metabolism, mTOR signaling pathway, and mismatch repair with FDR <0.25 ( Figure 3D). We further identified 25 immune related genes from Immport database for subsequent analysis ( Figure 3E, Table 2).

Identification of Hub TMB-Related Immune Genes and Associations of CNVs With Immune Infiltrates in SKCM
As the workflow in Figure S4, batch survival analysis indicated that four prognostic hub immune genes that highly associated with survival outcomes. Higher expression levels of CNTFR (ciliary neurotrophic factor receptor), CRABP2 (cellular retinoic acid binding protein 2), GAL (galanin and GMAP prepropeptide), and PAEP (progestogen-associated endometrial protein) were correlated positively with poor survival outcomes ( Figures 4A-D). Meanwhile, the prognostic value of four hub immune genes were validated in uveal melanoma cohort from TCGA. The results showed that higher expression of CNTFR, CRABP2, and PAEP were correlated positively with poor survival outcomes ( Figure S5). Moreover, we also analyzed the underlying relationship of the CNVs of four hub genes with immune infiltrates in SKCM microenvironment. Compared with the immune infiltration levels in samples with normal copy number of the signature, diverse forms of CNVs carried by four hub genes could commonly inhibit immune infiltrates, including CD8+ T cell, CD4+ T Cell, neutrophil cell, dendritic cell, macrophage, and B cell (Figures 5A-D).

Differential Abundance of Immune Cells in Two TMB Groups
Based on the newly developed CIBERSORT method, we calculated the proportions of 22 immune cells in each SKCM sample (Table S5) and exhibited the result in box plot, in which various colors represented different cell subsets ( Figure 6A). Besides, Wilcoxon rank-sum test revealed that infiltration levels   of regulatory T cell and dendritic activated cell in high-TMB group were lower than that in low-TMB group, while M1 and M2 macrophages showed the opposite trend ( Figure 6B).

Low Neutrophil and Macrophage Infiltrates Confer Poor Survival Outcomes
To investigate the underlying prognosis of immune cells, we performed the univariate analysis of infiltration levels of six immune cells associated with OS. The results indicated that lower infiltration levels of B cell, CD8+ T cell, neutrophil, and dendritic cell correlated with poor survival outcomes in SKCM ( Figure 7A). Furthermore, we conducted the multivariate Cox regression model and the data showed that lower infiltration levels of macrophage and neutrophil were risk factors for SKCM ( Table 3).

Construction and Assessment of TMBPI for SKCM
Given the relationship between alteration of immune signature with lower immune infiltrates and poor prognosis, we further evaluated the predictive accuracy of the four hub TMB-related immune genes. Based on the multivariate Cox regression analysis, we constructed the TMBPI as the following formula: TMBPI = 0.011581 × CNTFR + 0.006230 × GAL + 0.002504 × CRABP2 + 0.000261 × PAEP. Then, we calculated the TMBPI for SKCM patients and divided them into two TMBPI levels using the median value as the cutoff ( Table S6). The ROC curve of 3year OS prediction was drawn to assess the predictive accuracy with AUC = 0.664 ( Figure 7B). Meanwhile, Kaplan-Meier analysis indicated that SKCM patients with high TMBPI revealed worse survival outcomes compared with low TMBPI, which warranted further investigation and larger samples to validate ( Figure 7C).

DISCUSSION
In recent years, there has been a rapid development of immunotherapy for the treatment of advanced melanoma. Nowadays, PD-1 monotherapy is the standard first line therapy for advanced cutaneous melanoma, with efficacy, toxicity, and their correlations well established (30). However,  the data from real world showed that durable responses and favorable long-term outcomes are limited to a fraction of patients (31). Thus, many studies are designed to find predictive biomarkers for immune responses. Indeed, biomarkers remain the major challenge for immunotherapy, not only to identify patients who could respond to treatment but also to avoid unnecessary costs and serve toxicities for non-responders, including hyperprogression under ICB therapy, and to identify candidate patients for combination treatment. TMB and TILs, the novel biomarkers to predict immune responses, had been demonstrated their efficacy in several types of cancer, such as lung cancer, colorectal cancer, cutaneous melanoma, and so on (15,17,22,25,32). Nevertheless, few relevant researches had focused on the association of TMB with immune infiltrates and their prognostic role in SKCM.
In our study, we investigated the status of TMB in SKCM. The evidence from the analysis of the landscape of mutation profiles in our cohort indicated that 93.36% of patients contain diverse types of mutation. This was consistent to a previous study in which cutaneous melanoma have high mutation load and a predominant C>T nucleotide transition signature attributable to ultraviolet radiation (33). Analysis of correlations between TMB and clinicopathological characteristics showed that higher TMB level correlated with higher age, lower AJCC-T stage, lower AJCC-N stage, and early pathological stage. These data were in accordance with similar results in other clinical trials (NCT01295827 and NCT01866319) that patients with higher age tend to be more sensitive to ICB therapy (34,35), indicating that TMB is a promising biomarker for response to ICB therapy in SKCM to some extent. Besides, SKCM patients with early pathological stage may be more sensitive to ICB therapy. Inconsistent with previous investigations (26), no significant difference of TMB level was observed with AJCC-M stage in this study. This disparity may be mainly due to unbalanced stages. In our cohort, the frequency of samples with M1 stage was less than 10%.
We also analyzed the relationship between TMB and immune infiltrates in SKCM. We found that most of immune-related genes were downregulated in the lower-TMB group. Meanwhile, the Treg cells were inclined to be upregulated in the lower-TMB group, indicating that high TMB may inhibit immune cell infiltration in the tumor immune microenvironment. Moreover, CNVs of the four hub immune genes were related to immune cell infiltration in SKCM microenvironment. Immune infiltrates including CD8+ T cell, CD4+ T cell, neutrophil cell, dendritic cell, macrophage, and B cell could be inhibited by diverse forms of CNVs carried by hub genes. It implies that TMB have a significant impact on immune cell infiltration, resulting in alteration of tumor immune microenvironment in SKCM. When compared survival prognosis between the higher-TMB group and lower-TMB groups, we found that the higher-TMB group had longer survival time than the lower-TMB group. These data indicated that TMB is correlated to survival prognosis in SKCM, and that the reason leading to this correlation could attribute to the marked differences in immune infiltrates.
In contrast, other researches discovered that the patients in higher-TMB levels are inclined to have poor survival outcomes compared with lower-TMB levels in some types of cancer, including colorectal cancer, clear cell renal cell carcinoma, head and neck squamous carcinoma (26,29). Seemingly, the association between TMB and survival prognosis exhibited discrepancy in different types of cancer. The disparity may be mainly due to the most of SKCM patients were likely treated with immunotherapy. The data from our cohort and other cancer cohorts implicated that higher-TMB patients could gain a better prognosis than that in lower-TMB patients if treated with immunotherapy, otherwise higher-TMB patients would have poor prognosis compared with lower-TMB patients.
Finally, a prognostic model (TMBPI) was constructed using four hub immune genes which can be used for predicting survival outcomes in SKCM. Patients with high TMBPI revealed poor survival outcomes compared with that with low TMBPI. However, the AUC of the ROC curve was only 0.664 and further large-scale researches are required for verification and modification before clinical application.
There are some limitations, unresolved concerns, and potential perspectives in our study. First, we analyzed the association of TMB with immune infiltrates and their prognostic value using the data from TCGA cohort. The testing cohort from our clinical practice were lacking to evaluate the prognostic role of TMB and its relationship with immune infiltrates. Second, we used CIBERSORT technique to calculate the proportions of immune cells in every patient, which may simply analyze the cell composition on a large scale.  However, further flow cytometry and immunohistochemistry are necessary to validate the results. Besides, the genes with "Multi_hit" were not further assessed to figure out which of those genes were the actional mutations. Finally, the basic experiments were lacking to validate the association between four immune genes signature and immune infiltrates.
In summary, we revealed that higher TMB levels correlated with better clinical outcomes and may inhibit immune cell infiltrates in SKCM. Our study also identified four hub TMBrelated immune genes, and the CNVs of four hub genes conferred lower immune cells infiltrates. In addition, the prognostic model we constructed indicated that higher TMBPI conferred poor survival outcomes, which warranted further investigation and larger samples to validate.

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

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Medical Ethics Committee of Henan Provincial People's Hospital. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.