Immune cytolytic activity in skin cutaneous melanoma
We assessed the intratumoral immune cytolytic activity (CYT) in melanoma patients, measuring the TPM values of their GZMA and PRF1 . Table S1 lists the exact TCGA-SKCM dataset samples that were included in our analysis, along with their clinical information. We then stratified patients by defining skin melanomas in upper quartile of the cytolytic index, as CYT-high, and those in the lower quartile, as CYT-low. GZMA and PRF1 were tightly co-expressed across both primary and metastatic tumors (Spearman rank correlation, rho ~ 0.9) (Fig. 1a).
We and others have previously shown that CYT is markedly higher in skin melanoma relative to the normal skin [14, 36]. Here, we found that metastatic melanomas showed significantly higher cytolytic activity compared to primary tumors (p = 8.11e-06). Importantly, this was not observed between primary and metastatic tumors falling within the same extreme CYT percentiles, suggesting the existence of two main subgroups of tumors: CYT-high and CYT-low (Fig. 1b).
At the protein level, it was evident that neither toxin was highly expressed in skin tumors. In specific, GZMA protein expression (nuclear or cytoplasmic/membranous) was medium only in one skin melanoma; whereas, it was low in six and absent in five out of twelve skin tumors. The corresponding intensity scores were as follows: moderate, 2/12; weak, 7/12; and negative, 3/12. The absence of PRF1 protein was more evident, as it was not detected in any of the 11 skin melanomas (Fig. 1c and Table S2).
As expected, metastatic melanomas had significantly higher total number of mutations against primary melanomas (p = 0.00178) (Fig. 1d). However, when we assessed the CYT-high and CYT-low skin tumors separately, the metastatic subgroup did not exhibit any difference in the mutation load (p = 0.772) and this was not correlated with the cytolytic index. Among primary tumors on the other hand, CYT-high melanomas accumulated significantly more mutations (p = 0.0186), and the mutational burden was significantly correlated with the tumors’ cytolytic activity, suggesting a primary association between cytolytic activity and mutation load in these tumors (Fig. 1e-f).
In tune with the mutagenic role of UV radiation in melanoma, the most prevalent mutational signatures that we detected had higher cosine-similarity with COSMIC's signatures 7 (CC > TT dinucleotide mutations at dipyrimidines due to UV exposure) and 11 (strong transcriptional strand-bias for C > T substitutions, resembling treatment with alkylating agents). This preference for signature 7 (RSS < 1e-03, cosine similarity > 0.99) did not differ between the two cytolytic subgroups, whatsoever (Figure S1a-b). They also resembled with signature 30 (C > T, unknown aetiology, observed in a small subset of breast cancers) and to a less degree (cosine similarity, 0.724) to mutational signature 6 (defective mismatch repair system) (Fig. 1g-h).
Kataegis is equally distributed across different cytolytic subgroups of skin melanoma
We hypothesized that mutation showers (or kataegis) are associated with cytolytic-high cutaneous melanomas. Kataegis were recently uncovered by whole-genome sequencing of B cell and non-hematopoietic tumors [58, 60, 63–65]. We identified 74 kataegic sites across 26 tumors, associated with 1,567 mutations. Five (19.2%) of these tumors were CYT-high metastatic SKCMs and six (23.07%) CYT-low (3 metastatic and 3 primary tumors), suggesting that kataegis is equally distributed across skin tumors of different cytolytic activity.
Importantly, 1,395 of the mutations (89%) were C > T transitions, which is consistent with the notion that kataegis results from DNA replication over cytidine deamination of resected DNA [65, 66]. Considering the established role of the APOBEC cytidine deaminases in kataegis , these data support a partial APOBEC involvement in cutaneous melanomas, irrespective of their immune cytolytic index.
We next investigated whether CYT associates with transcriptional or genomic differences between different cytolytic subgroups in primary and metastatic skin melanomas.
Cyt Correlates With Distinct Gene Sets In Melanoma
Initially, we mapped the reads to the GRCh37 (hg19) human reference genome and summarized them the gene level using Rsubread  producing a matrix of counts. The read counts were converted to log2-counts-per-million (logCPM) and voom transformed using limma  to detect the significantly differentially expressed genes in skin melanoma (primary and metastatic) vs the normal skin tissue (Figure S2). Using principal-components analysis (PCA), it was evident that the differentially expressed genes in metastatic skin melanomas could discriminate them better from the normal skin, compared to the differential genes in primary tumors. We detected 1,054 (11.5%) co-deregulated genes (adj.p < 0.001) between primary and metastatic skin melanomas, but also 13 (0.1%) genes that were deregulated only in primary and 8,121 (88.4%) genes in metastatic skin melanomas (Table S3). Importantly, PRF1, GZMA, NKG7, SLA2, GBP5 and CD2 were among the top 20 upregulated genes both in primary and metastatic skin melanomas. In primary SKCM the top 20 deregulated genes also included CD8A, LAG3, PDCD1, ITGAL, TRAC, TIGIT, TNFRSF9, CD247, CD27, SIT1, CD8B, CCR5, SIRPG and CD3D (Table S3 and Figure S3). On the other hand, the top 20 upregulated genes among metastatic tumors also included IRF1, CD74, FASLG, JAKMIP1, C1QA, IFNG, TNIP3, C1QB, APOL3, C1QC, CCL5, TRGV10, CRTAM and CD8A (Table S3 and Figure S3).
To better evaluate the sources of gene expression across all SKCM tumors, we calculated tumor purity and ploidy, using ABSOLUTE . Overall, CYT-low tumors had higher purity compared to the CYT-high ones; whereas, the ploidy values of most tumors clustered around 2.4–2.7 (genome-wide duplication), without differences between the two cytolytic subsets (Figure S4).
Using GSVA on RNA-seq data extracted from the TCGA-SKCM dataset and the “C2: curated gene sets” from MSigDB v6.1 (containing 4,738 gene sets) , we assessed enrichment within the CYT-high primary skin melanomas, for the immune-related gene sets BIOCARTA_TCRA_PATHWAY, REACTOME_TRANSLOCATION_OF_ZAP_70_TO_IMMUNOLOGICAL_SYNAPSE, REACTOME_PD1_SIGNALING, BIOCARTA_TCAPOPTOSIS_PATHWAY, CHAN_INTERFERON_PRODUCING_DENDRITIC_CELL and REACTOME_ENDOSOMAL_VACUOLAR_PATHWAY. On the other hand, CYT-low primary tumors were statistically enriched for DNA replication and DNA-repair gene sets, including REACTOME_REPAIR_SYNTHESIS_FOR_GAP_FILLING_BY_DNA_POL_IN_TC_NER, REACTOME_POL_SWITCHING, REACTOME_LAGGING_STRAND_SYNTHESIS, REACTOME_PROCESSIVE_SYNTHESIS_ON_THE_LAGGING_STRAND, REACTOME_DNA_STRAND_ELONGATION, REACTOME_REMOVAL_OF_THE_FLAP_INTERMEDIATE_FROM_THE_C_STRAN, REACTOME_UNWINDING_OF_DNA, BIOCARTA_MCM_PATHWAY, REACTOME_ACTIVATION_OF_THE_PRE_REPLICATIVE_COMPLEX, KALMA_E2F1_TARGETS, CROSBY_E2F4_TARGETS, among others (Fig. 2).
CYT-high metastatic melanomas were also enriched for immune-related gene sets, including KEGG_AUTOIMMUNE_THYROID_DISEASE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_ALLOGRAFT_REJECTION, BIOCARTA_CTL_PATHWAY, BIOCARTA_TCYTOTOXIC_PATHWAY, BIOCARTA_THELPER_PATHWAY, BIOCARTA_THELPER_PATHWAY, BIOCARTA_TCAPOPTOSIS_PATHWAY, BIOCARTA_TCRA_PATHWAY, REACTOME_PD1_SIGNALING, and BIOCARTA_DC_PATHWAY (Fig. 2); whereas, CYT-low metastatic tumors were enriched for non-immune-related gene sets (TESAR_ALK_TARGETS_HUMAN, MYLLYKANGAS_AMPLIFICATION, R_KINOME_KUMAMOTO_RESPONSE_TO_NUTLIN_MONTERO_THYROID_CANCER_POOR_SEMBA_FHIT_TARGETS_DN, CROSBY_E2F4_TARGETS
REACTOME_UNWINDING_OF_DNA, ZERBINI_RESPONSE_TO_SULINDAC_OHASHI_AURKA_TARGETS, etc.).
Interferon-gamma (IFN-γ) is an important cytokine produced by activated T cells, NK and NK T cells, in the tumor microenvironment, and it plays a key role in orchestrating innate and adaptive antitumor immune response . To gain further insight into the relationship of IFN-γ and the cytolytic index in skin melanoma, we explored two IFN-γ signatures that were recently proposed to predict patient response to Pembrolizumab . These signatures contain IFN-γ–responsive genes related to antigen presentation, chemokine expression, cytotoxic activity, and adaptive immune resistance (“IFN-γ” signature: IFNγ, IDO1, CXCL9, CXCL10, HLA-DRA, STAT1; and “expanded immune” signature: CD30, IDO1, CIITA, CD3E, CCL5, GZMK, CD2, HLA-DRA, CXCL13, IL2RG, NKG7, HLA-E, CXCR6, LAG3, TAGAP, CXCL10, STAT1, GZMB). Both gene signatures were markedly overexpressed in CYT-high tumors, corroborating the enrichment that we found for the similar gene set CHAN_INTERFERON_PRODUCING_DENDRITIC_CELL across these tumors (Figure S5), and the high discriminatory value of IFN-γ–related gene signatures, in enriching response rates to Pembrolizumab .
Taken together, the above results provide further proof that the tumor microenvironment in CYT-high skin melanomas is more inflamed and immunogenic compared to that in CYT-low tumors.
Other than GZMA and PRF1, the immune checkpoints CTLA-4, PDCD1 (PD1), CD274 (PD-L1), PDCD1LG2 (PD-L2) and IDO1 were significantly overexpressed in CYT-high (primary and metastatic) SKCM tumors, underlying the immunosuppressive microenvironment of these tumors (Fig. 2c-d).
To assess the relationship between CYT markers and immune checkpoint molecules, we run correlation analysis of their gene expression levels. Both cytolytic genes were notably correlated with the expression of at least five different inhibitory checkpoints in SKCM tumors (p < 0.0001), corroborating that combination therapies with immune checkpoint inhibition should effectively be used to overcome resistance and broaden the clinical benefit for these patients (Fig. 2e).
In addition, we examined the expression of the immunoinhibitor LAG3, the immunostimulators CD70, lipoteichoic acid (LTA) , ecto-5′-nucleotidase (NT5E)  and ectonucleoside triphosphate diphosphohydrolase 1 (ENTPD1) , as well as several regulatory cytokines and chemokines, known for their pro- and anti-inflammatory roles within the tumor microenvironment [14, 35]. The expression of these molecules was also compared between the two cytolytic subgroups in primary and metastatic SKCM. In specific, we focused on markers for activated CD8 + T cells (NKG7, CD3E, GZMA, GZMH, GZMK), MDSCs (CD2), activated dendritic cells (UBD, C1QB, C1QC), NK cells (FASLG and FAS), and interferon-stimulated chemokines that attract T cells (CXCL9, CXCL10, CXCL11 and CXCL13) . All these genes were significantly upregulated in CYT-high (primary and metastatic) tumors, supporting the existence of an inflamed microenvironment in them (Figure S6).
Cell type fraction analysis revealed that CYT-high (primary and metastatic) tumors are significantly enriched in B cells, M1 macrophages and CD8 + T cells; whereas CYT-low tumors contain significantly higher levels of monocytes, NK cells and CD4 + T cells (Fig. 2d). Further CIBERSORT analysis revealed that the majority of (primary and metastatic) CYT-high tumors were significantly enriched in γδT cells, follicular helper T cells (Tfh) and Tregs, compared to their CYT-low counterparts (Figure S7). These cells were previously shown to participate in cancer immunosurveillance , autoimmunity  and immunosuppression .
Cyt Correlates With Discrete Mutational Events In Melanoma
We next focused on the exome-seq data of the TCGA-SKCM and detected the significantly mutated genes in each cytolytic subgroup of skin melanoma. CYT-low primary tumors were significantly associated with mutations in BRAF, DMRT3, GSTA5, TP53, CRYBA4, STRA13, PRAMEF12, RNF32, SLC1A6, TEKT2 and NRAS among others; whereas, CYT-high primary tumors with a totally different group of genes, including LCK, GSTSF1L, BMF, CSTL1, DRGX, ENTPD3, CAMK4, GPR151, LTF and others (Fig. 3a). This discrepancy hints that the two cytolytic subsets are correlated with different somatic mutations; but we should also take into consideration the small sample number of CYT-high primary tumors (n = 10) compared to the CYT-low ones (n = 31).
In contrast, the SMGs in both cytolytic subgroups of metastatic melanomas included previously described  driver mutations in oncogenes and tumor suppressors (NRAS, BRAF, CDKN2A, TP53 and PTEN), as well as COL4A4.
CYT-low metastatic melanomas were associated with non-silent mutations in GALNTL5, RGPD4, DNAJC5B, DEFA3, IRF2, COL2A1, VEGFC, WFDC11, PROL1, CD80, ADAM18, APOBEC3H and CSHL1; whereas CYT-high metastatic skin tumors were significantly mutated in PPP6C, DSG1, DEFB112, SPATA16, STARD6, RARRES2 PDE1A, IQCF3, KLK8, CDKN2A, SPANXN5, C16orf90, CD48, TXNDC3 and RAC1 (Fig. 3a).
The fraction of C > T transitions at dipyrimidines was the highest across all skin melanomas. The frequency of specific substitutions did not differ between CYT-high and -low metastatic tumors. However, across the cohort of primary tumors, CYT-high melanomas had more C > T transitions against CYT-low tumors (p < 0.01); but this could probably be due to the difference in sample number (103 primary vs 368 metastatic skin melanomas) (Fig. 3b).
Likewise, there was no association between CYT and BRAF, NRAS, TP53 or NF1 mutations, indicating that immune-related alterations within the tumor microenvironment, and therefore immunotherapies, are independent of the tumor’s genotype (Figure S8).
Cyt Associates With Different Structural Changes In Melanoma
Skin melanoma is characterized by increased chromosomal instability (CIN) with extensive gains and losses, which associate with poor patient prognosis [81–83]. However, whether these chromosomal aberrations correlate with the cytolytic index, is unknown. Therefore, we assessed the somatic copy number aberrations (SCNA) within each cytolytic subtype in primary and metastatic skin melanomas.
CYT-low primary tumors had recurrent amplification in locus 1p12 (NOTCH2) and deletions in loci 9p22.4 (JAK2, PDL1/CD274, PDL2/PDCD1LG2), 10q23.1 (miR346), 11q25 (ETS1), 15q15.1 (RAD51), 16q23.3 (CDH13) and 6q27 (CCR6, RNASET2, FGFR1OP). On the other hand CYT-low metastatic melanomas had recurrent amplifications in loci 5p15.33 (TERT), 6p25.1 (CDYL), 7p22.2 (RAC1), 12q15 (MDM1) and recurrent deletions in 5q31.2 (CTNNA1, FGF1), 11q22.3 (FDX1, RDX, ZC3H12C), and 15q14 (RASGRP1, CSNK1A1P1, SPRED1).
CYT-high primary melanomas on the other hand, did not have recurrent copy number amplifications above the threshold (Q value = 0.25), but rather had loses at 9p21.3 (CDKN2A/B), 9q34.11 (CDK9) and 16q24.2 (IRF8). CYT-high metastatic melanomas had recurrent copy number amplifications in loci 1p12 (NOTCH2), 1q44 (OR2M4), 7q36.1 (KCNH2), 11q13.3 (FGF3/4), 12p11.23 (MED21) and 12q13.3 (R3HDM2) and deletions in 1p22.1 (RHOC, NRAS), 9p21.3 (CDKN2B), 10q23.2 (miR346), 11q23.3 (ATM, CHEK1, ETS1) and 15q15.3 (CASC4).
CDKN2A/B deletions (9p21.3) were common between CYT-high and CYT-low melanomas (both primary and metastatic). Also, amplifications in loci 1p12 (NOTCH2), 3p13 (MITF), 13q13.3 (FGF3/4, CTTN) and 22q13 (MLK1) were found in both cytolytic subtypes of metastatic melanoma (Fig. 4a). Overall, metastatic melanomas had significantly more recurrent SCNAs compared to primary tumors (602 amplifications and 3,920 deletions in metastatic SKCM relative to 99 amplifications and 1,632 deletions in primary SKCM), but also CYT-low tumors had significantly more SCNAs compared their CYT-high counterparts (Fig. 4b).
Taken together, the above findings show that specific types of somatic mutations and copy number aberrations are characteristic for each cytolytic subgroup in primary and metastatic melanoma.
To eliminate the possibility of assessing lower confidence of detecting somatic mutations and SCNAs due to tumor cellularity, we performed ABSOLUTE analysis  and found no variance in the calculated cellularity estimates between CYT-high and –low SKCM samples. MATH scores ranged from 14.70-47.17 in primary and 14.59–62.22 in metastatic SKCMs. Neither the total mutation load, nor the total copy number events, correlated with tumor heterogeneity in melanoma. In addition, the MATH scores were similar between the two cytolytic subgroups (in primary SKCM, CYT-high vs. CYT-low, 31.04 ± 11.52 vs 29.86 ± 12.97, p > 0.05; in metastatic SKCM, 34.43 ± 12.02 vs 29.94 ± 13.31, p > 0.05), signifying that the dissimilarities in copy number and mutational load, are not due to a variable intra-tumor heterogeneity (Fig. 4c). Therefore, we propose that distinctive mutational and structural changes can discriminate skin melanomas of different cytolytic activity.
Overexpression of GZMA and PRF1 synergistically affects patient overall survival
To determine whether CYT-high is a good prognostic indicator, we explored the overall survival of skin melanoma patients having high or low expression levels of GZMA and PRF1, using SynTarget .
High expression of each cytolytic gene, individually, was positively correlated with the overall survival in melanoma patients (p = 0.02). The subgroup of non-metastatic melanoma patients having both genes over-expressed showed a significant positive effect in survival compared to the remaining patients (“other”) (p = 0.0282); whereas, simultaneous low levels of both genes shifted significantly (p = 0.0048) towards a negative effect (Fig. 4d). This observation implies that the overexpression of both GZMA and PRF1 genes, can synergistically affect the overall survival of skin melanoma patients.
Chromothriptic Events In Cytolytic Subsets Of Skin Melanoma
Local chromosome shattering (chromothripsis) is a mechanism proposed to cause clustered chromosomal rearrangements, following chromosomal breaks at multiple locations. Chromothripsis has been detected in 2–3% of cancers  and involves impaired DNA repair . Such genomic rearrangements may drive the development of cancer through the deletion of tumor suppressor genes or an increase in copy number of oncogenes, among other mechanisms. The occurrence of chromothriptic events in skin melanoma, and their effect on genes associated with checkpoint inhibition and immune cytolytic activity is still unclear.
Overall, we detected chromothriptic events of various sizes, in 211/470 (44.89%) skin melanomas (Fig. 5a, d-f). Among these, 25 belonged to the CYT-high and 41 to the CYT-low subgroups (average CNA Status change times, 50.6 ± 35.17, CYT-high vs 41.12 ± 21.59, CYT-low SKCM). Notably, we observed distinct patterns of chromothriptic events between the two cytolytic subgroups, the majority of which were harbored mainly in CYT-low skin melanomas (41/66, 62% of the chromothriptic events, CN change times ≥ 20 and log10 of LR ≥ 8) (Fig. 5b) and affected more chromosomes in CYT-low tumors (Fig. 5c). Interestingly, the chromosomal regions among CYT-low skin melanomas that harbored chromothripsis, contained a higher number of cancer genes, including KRAS, NOTCH2, BCL9, CCND1 (gains) and BRAF, NRAS, PAX3, ATM and CD274 (losses) (306 gains and 209 losses in total, among CYT-low tumors vs. 156 gains and 75 losses in total, among CYT-high tumors (Fig. 5b and Table S4).
The existence of different chromothriptic patterns has been proposed to be due to different mechanisms, including dysfunction of telomeres and DNA damage in micronuclei [54, 86–89].
The telomere regions were not affected by chromothripsis in the vast majority of events (93,521 out of 93,535 events, 99.9%) across all skin melanomas. Focusing on the two cytolytic subgroups though, the telomere region affected one CYT-high and seven CYT-low skin tumors out of the 14 events, in total (Fig. 5b).
Chromothriptic events not directly affecting telomere regions can also result from telomere fusions, since the genomic regions included in chromatin bridges can be distant from the telomeres, depending on the structure of the dicentric chromosomes formed in telomere crisis [54, 86]. In more than one fourth of the tumors (58/211; 27.5%), the centromere was included in the segment that was affected by chromothripsis, including 16 CYT-high and 68 CYT-low tumors (Fig. 5b).
Overall, our findings support the existence of a broader pattern of chromothripsis across cytolytic low skin melanomas.
Increased neoantigen load associates with high cytolytic levels in primary skin melanomas
Cancer neoantigents arise from mutated peptides and favorably drive T-cell recognition of cancer cells [90, 91]. They are therefore, an attractive immune target because their selective expression on cancer cells can minimize immune tolerance .
It was of interest to confirm whether cancer neoantigents associate with high cytolytic levels in skin melanoma. Analyzing mutation data from the Cancer Immunome Database (TCIA) , it was evident that primary CYT-high skin tumors contained higher number of mutations and neoantigens compared to their CYT-low counterparts. This was not evident, however, in metastatic CYT-high tumors (Fig. 6a).
To investigate this further, we predicted the missense mutations that could potentially function as T-cell neoepitopes in skin melanomas, using antigen.garnish. Our results confirmed that CYT-high primary (and not metastatic) tumors are significantly enriched in (classically and alternatively defined) neoantigens, but these, were not correlated with a high cytolytic activity (Fig. 6b). The tumors’ MATH scores also did not correlate with their neoantigen load. Overall, these data support that a high cytolytic activity is driven by a high mutational/neoantigenic load only in primary skin melanomas.
Prediction of skin melanoma patients’ response to immune checkpoint inhibition
Since just a subgroup of skin melanoma patients responds to immune checkpoint inhibition, the need to elucidate the mechanisms of resistance and predict markers, is high. TILs, the expression of PD-1 or PD-L1, the mutational load , or clonal neoantigens , have all been proposed as markers; however, none of them has been fully validated, yet .
Hypothesizing that CYT-high skin melanomas have higher immunophenoscore due to increased immunogenicity, resulting in better prognosis and response to therapy, we analyzed two data sets containing skin melanoma patients treated with anti-CTLA-4  and anti-PD-1 inhibitory molecules , and used each patient's IPS to predict their response.
Assessed globally, CYT-high skin melanomas (both primary and metastatic) had significantly higher HLA levels compared to the intratumoral mean expression. In contrast, CYT-low melanomas were characterized by downregulation of several immune checkpoints, compared to the intratumoral mean expression. Importantly, we could observe similar patterns across patients treated with anti-CTLA-4 alone or combined with anti-PD-1/PD-L1/2, or anti-PD-1/PD-L1/2 monoclonal antibodies, alone (Fig. 7a-b). In addition, CYT-high SKCMs contained significantly higher numbers of cytotoxic cells (CD8 + T cells, γδT cells, NK cells), and lower numbers of MDSC and Treg cells.
As expected, CYT-high SKCM patients had a markedly higher IPS (and consequently, an expected clinical benefit) compared to CYT-low patients who did not receive checkpoint inhibition. Interestingly, the IPS scores were significantly higher across all CYT-high tumors (both primary and metastatic), upon treatment with CTLA-4 or PD-1/PD-L1/2 blockers, or a combination of both (p < 0.0001) (Fig. 7c).
Our data indicate that the IPS has a predictive value in CYT-high melanoma patients who received CTLA-4 and PD-1/PD-L1/2 inhibition therapy, and are in accordance with previous observations that patients with higher levels of tumor cytolytic activity, and expression of immune checkpoints, benefited more from the corresponding immune checkpoint blockade .