This is a preprint. Preprints are preliminary reports that have not undergone peer review. They should not be considered conclusive, used to inform clinical practice, or referenced by the media as validated information.
Research
Distinct Genomic Features Across Cytolytic Subgroups in Skin Melanoma
Apostolos Zaravinos
Department of Basic Medical Sciences, College of Medicine, Member of QU Health, Qatar University, 2713 Doha, Qatar.
Corresponding Author
ORCiD: https://orcid.org/0000-0003-4625-5562
License:
This work is licensed under a Creative Commons Attribution 4.0 International License. Read Full License
Background: Skin melanoma is a highly immunogenic cancer with extensive genetic and transcriptional diversity. Despite the enthusiastic clinical results seen in advanced-stage metastatic melanoma patients treated with immune checkpoint inhibition, a subgroup of them later relapse and develop acquired resistance. The intratumoral immune cytolytic activity (CYT) reflects the ability of cytotoxic T cells and NK cells to eliminate cancer cells, and is associated with improved patient survival. Here, we questioned whether CYT associates with different genomic profiles in skin melanoma.
Methods: We explored the TCGA-SKCM dataset and stratified patients to distinct subgroups of cytolytic activity. The tumor immune contexture, somatic mutations and recurrent copy number aberrations were calculated using quanTIseq, MutSigCV and GISTIC2. Chromothriptic events were explored using CTLPScanner and cancer neoepitopes were predicted with antigen.garnish. Each tumor's immunophenoscore was calculated using Immunophenogram.
Results: Metastatic skin melanomas had significantly higher CYT compared to primary tumors. We also assessed enrichment for immune-related gene sets within CYT-high tumors; whereas, CYT-low tumors were enriched for non-immune related genes sets. In addition, distinct mutational and neoantigenic loads, primarily composed of C>T transitions, along with specific types of copy number aberrations, characterized each cytolytic subgroup. We found a broader pattern of chromothripsis across CYT-low tumors, where non-telomeric chromosomal regions harboring chromothripsis, contained a higher number of cancer genes. CYT-high patients had markedly higher immunophenoscore and should consequently, display an expected clinical benefit compared to CYT-low patients who either received or not, checkpoint inhibition.
Conclusions: Overall, our data highlight the existence of distinct genomic features across cytolytic subgroups in skin melanoma.
Keywords
skin melanoma, cytolytic activity, immune checkpoints, chromothripsis, kataegis, chromosomal instability, neoantigen load, immunophenoscore, cell type fractions
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Cutaneous melanoma is a very aggressive and highly immunogenic cancer with extensive genetic and transcriptional diversity [1–8]. The tumor usually results after exposure to UV radiation, and has the highest rate of somatic mutations and neoantigens among all cancer types [9]. The Cancer Genome Atlas (TGCA) Network recently classified skin melanoma into mutant BRAF, mutant RAS, mutant NF1, and Triple-WT genetic subtypes, based on the pattern of the top mutated genes [10].
The tumor’s microenvironment comprises a very heterogeneous cell population, including fibroblasts, lymphocytes, macrophages and other immune cells, as well as, adipocytes and cells that form the structural elements of skin blood vessels dashed in the extracellular matrix [11]. A high content of tumor-infiltrating lymphocytes (TILs) in the tumor’s microenvironment associates with a favorable prognosis and improved overall survival of the patients [12, 13]. An increased immune cytolytic activity (CYT), defined by the attempt of the cytotoxic T cells and natural killer (NK) cells to eliminate cancer cells via the secretion of GZMA and PRF1, has also been associated with improved patient survival [14].
With the advent of immunotherapies, the tumor’s management has shifted from cytokine-based treatment to immune checkpoint inhibition, primarily of cytotoxic T-lymphocyte-associated antigen-4 (CTLA-4) and programmed cell-death protein 1 (PD-1) or its ligands (PD-L1 and PD-L2) [15–18]. Monoclonal antibodies including Ipilimumab, Tremelimumab (anti-CTLA-4), Pembrolizumab, Nivolumab, Cemiplimab (anti-PD-1), Atezolizumab, Avelumab and Durvalumab (anti-PD-L1), among others currently being on clinical trial [19], have enthusiastically provided long-lasting responses and improved survival in patients with advanced-stage metastatic melanoma [16].
Early clinical results have also been stated with IDO inhibitors [20] and Treg depletion targeting surface CD25 [21, 22]. Many other immunotherapeutic agents are currently under study, targeting a range of positive or negative immunoregulatory molecules within the tumor microenvironment. These, include agonistic antibodies against 4-1BB, Ox40, Inducible T-cell COStimulator (ICOS), and CD40; blocking antibodies against Lymphocyte-activation gene 3 (LAG3), B7-H3, B7-H4, Tim3, and killer inhibitory receptors (KIRs); as well as cytokines IL-7, IL-15 and IL-21 [19, 23–32], among others.
Due to the complexity of immune regulation in vivo, combination immunotherapies are expected to provide a better therapeutic benefit for the patients [19]. Nevertheless, a subgroup of responders to immunotherapy, later relapse and develop acquired resistance; whereas, others do not respond at all (primary resistance) [33]. Therefore, a better understanding of the underlying resistance mechanisms is emergently needed.
Here, we questioned whether CYT associates with different genomic profiles in skin melanoma. We investigated how the immune landscape in these tumors relates to somatic mutations, mutational signatures, copy number aberrations, chromothriptic events, the expression of immune checkpoints or other immune-related markers, and the presence of different types of immune cells within the tumor microenvironment of primary or metastatic tumors. We also predicted how skin melanoma patients of different cytolytic activity respond to immune checkpoint blockade therapy. Overall, we provide enough evidence of the existence of distinct genomic features across cytolytic subgroups of skin melanoma.
Skin melanoma (SKCM) data extraction
We extracted the clinical information, Mutation Annotation files (MAF) and mRNA-Seq ‘level 3’ data of a total of 470 primary and metastatic skin melanomas (103 primary tumors, 68 distant metastasis, 74 regional cutaneous or subcutaneous tissues and 222 regional lymph nodes), along with their matched peripheral blood, from the Cancer Genome Atlas’ (TCGA-SKCM dataset) GDC Data Portal (https://portal.gdc.cancer.gov/). Three samples were non-informative and were thus, excluded from subsequent analyses. In addition, we accessed each patient’s gene-level, zero-centered, focal copy-number data analyzed by GISTIC (v2.0.22) [34] from Broad GDAC Firehose (https://gdac.broadinstitute.org/).
Cytolytic Activity Calculation And Downstream Rna-seq Analysis
We calculated each SKCM patient’s levels of activity (CYT), as the geometric mean of GZMA and PRF1 [14, 35, 36]. GZMA leads to caspase-independent apoptosis, while PRF1 forms pores in the target tumor cells facilitating the entry of granzymes into them. Gene expression values were presented in Transcripts Per Million (TPM).
We then divided patients into two cohorts, the upper 25th quartile of the cytolytic index (CYT-high) and the lower 25th quartile (CYT-low), each with an identical admixture of histology-stage combinations. All subsequent comparisons were made between CYT-high and CYT-low (metastatic or primary) skin melanomas. P-values were adjusted using the Benjamini-Hochberg (BH) method.
We used gene set variation analysis (GSVA) to refine alterations in pathway activity [37] and clustered tumors hierarchically, with complete linkage as the distance metric. We compared gene expression between the two cytolytic subsets of melanomas, first by estimating the mean-variance relationship of the log-counts using voom [38] and then analyzing them with limma [39]. We considered differentially expressed genes, those having a BH-adjusted p-value < 0.1. Graphs were plotted with ggplot2.
Immunohistochemistry (IHC) data of GZMA and PRF1 protein expression were extracted from the Human Protein Atlas [40–42], and further analyzed. The antibodies used in IHC were as follows: rabbit pAb anti-GZMA, HPA054134, 1:200 dilution, Sigma-Aldrich, Atlas Antibodies Cat#HPA054134, RRID:AB_2682395; Antigen retrieval was performed using HIER pH6; and mouse mAb anti-PRF1, CAB002436, 1:10 dilution, Leica Biosystems, Cat#NCL-PERFORIN, RRID:AB_563955; Antigen retrieval was performed using HIER pH6.
Cell Type Fractions And Heterogeneity
We quantified the tumor immune contexture from the RNA-seq data, using quanTIseq (https://icbi.i-med.ac.at/software/quantiseq/doc/) [43]. Additionally, we used CIBERSORT to identify fractions of immune subpopulations in each cytolytic subgroup’s tissue (https://cibersort.stanford.edu/) [44]. Tumor purity and ploidy estimates were produced using ABSOLUTE [45] and plotted for each cytolytic subgroup of primary and metastatic tumors.
Detection Of Somatic Mutations And Copy Number Aberrations (scna)
We calculated the most significantly mutated genes (SMG) in each CYT subtype of primary and metastatic skin melanomas, using MutSigCV (v1.3.01) [46, 47] with an FDR (q-value) = 0.1 as threshold of significance. We further employed GISTIC (v2.0.22) [34, 48] to identify genomic regions with significant amplifications or deletions in each sample. A G-score considering the amplitude of the aberration, as well as the frequency of its occurrence across each tumor sample was assigned. We calculated the somatic copy number alterations in each sample by taking the sum of segment mean changes ≥ 0.6 and ≤ -0.4 between somatic and normal samples. Regions with FDR q-values < 0.1 were considered significant. Both mutational and copy number analyses were performed using the Broad Institute’s GenePattern platform (http://software.broadinstitute.org/cancer/software/genepattern).
Detection Of Chromothriptic Events
Chromothriptic events in CYT-high or low skin melanomas were investigated using CTLPScanner [49]. In brief, after downloading segmentation data (level 3) of single nucleotide polymorphism (SNP) arrays from the TCGA-SKCM dataset, we implemented DNA copy number segmentation data, using the circular binary segmentation algorithm [50], and the NCBI38/hg38 genome assembly. For the detection of chomothripsis or chomothripsis-like regions, we used copy number status change ≥ 20 times, log10 of likelihood ratio ≥ 8, minimum segment size = 10 kb, and signal distance between adjacent segments = 0.3. Shattered chromosomomal regions were visualized based on the signal value for genomic gains (≥ 0.15) or losses (≤ -0.15) and further highlighted. Finally, chromothripsis-located genes were annotated using the COSMIC database [51].
Consensus centromeric regions were extracted from http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/centromeres.txt.gz. Telomeric regions were defined as 20 kb from each side of the chromosomal ends. Intersection between chromothripsis events and centromeric or telomeric regions was defined as any overlap between the event coordinates and the region coordinates.
Histograms of copy number segment switches were drawn to identify events of chromothripsis that arise in conjunction with telomere crisis, i.e., a period of genome instability during tumorigenesis when depletion of the telomere reserve generates unstable dicentric chromosomes [52–54]. The patterns and distribution of copy number state switches indicate the emergence and temporal order of the major genomic rearrangement events.
Detection Of Cancer Neoepitopes And Immunophenoscores
We predicted the cancer neoepitopes using “antigen.garnish”, with the following MHC molecules, as previously described [35, 55, 56]: H-2-IAb, H-2-IAd, HLA-DPA1*01:03, HLA-DPA1*02:01, HLA-DPA1*02:01, HLA-DPA1*03:01, HLA-DPB1*03:01, HLA-DQA1*01:01, HLA-DQA1*01:02, HLA-DQA1*03:01, HLA-DQA1*04:01, HLA-DQA1*05:01, HLA-DQA1*05:01, HLA-DRB1*01:01. We classified the peptides that were predicted to bind MHC with high affinity (IC50 < 50 nM) or with greatly improved affinity compared to their wild-type counterparts, either as classically (CDNs) or alternatively defined neoepitopes (ADNs).
The Cancer Immunome Database (TCIA) (https://tcia.at/) was also queried for cellular composition of neoantigens among the two cytolytic subsets of skin melanomas [57]. The immunophenoscore (IPS) of each cytolytic subgroup’s patients was calculated using Immunophenogram, as previously described [57], and immunophenograms were constructed to visualize the immunophenotypes of each tumor sample. The IPS ranged from 0 to 10, according to the expression of MHC molecules, immunomodulators, effector cells (activated CD8 + T cells and CD4 + T cells, effector memory (Tem) CD8 + and CD4 + cells), and suppressor cells, including regulatory T cells (Tregs) and myeloid-derived suppressor cells (MDSCs) [57].
Mutational Signatures And Rainfall Plots
We extracted mutational signatures as previously described [58, 59] and compared them against the validated mutational signatures provided by COSMIC (http://cancer.sanger.ac.uk/cosmic/signatures). We further applied cosine similarity against COSMIC's validated signatures to identify the best matches within signatures. We reconstructed the mutational spectrum of CYT-high and CYT-low skin melanomas, using the residual sum of squares (RSS) to measure the efficiency of the reconstruction of the original mutational profile.
The genomic regions of each cytolytic subgroup that exhibited rainfall plots, i.e., contained ≥ 6 consecutive mutations with an average inter-mutation distance of ≤ 1 bp, were defined as kataegis [58, 60, 61].
Tumor Heterogeneity
Tumor heterogeneity in SKCM tumors was inferred by clustering variant allele frequencies (VAF). The extent of each skin melanoma’s intra-tumor heterogeneity was quantitatively calculated by the width of the VAF distribution. We assigned a mutant-allele tumor heterogeneity (MATH) score to each SKCM sample, as previously described [35].
Patient Survival And Synergistic Target Analysis
The overall survival of the melanoma patients was performed using data extracted from the human skin cutaneous melanoma dataset (TCGA-SKCM). Kaplan–Meier curves using log rank (Mantel Cox) test estimated the patient survival. We examined whether GZMA and PRF1 act synergistically on patient survival, using SynTarget [62].
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 [14]. 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 [67], 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 [68] producing a matrix of counts. The read counts were converted to log2-counts-per-million (logCPM) and voom transformed using limma [69] 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 [45]. 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) [70], 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 [71]. 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 [72]. 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 [72].
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) [73], ecto-5′-nucleotidase (NT5E) [74] and ectonucleoside triphosphate diphosphohydrolase 1 (ENTPD1) [75], 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) [76]. 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 [77], autoimmunity [78] and immunosuppression [79].
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 [80] 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 [45] 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 [62].
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 [84] and involves impaired DNA repair [85]. 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 [92].
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) [57], 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 [93], or clonal neoantigens [91], have all been proposed as markers; however, none of them has been fully validated, yet [94].
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 [95] and anti-PD-1 inhibitory molecules [96], 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 [95].
A better understanding of the immunobiology of cancer immunosurveillance and immunoediting in skin melanoma, will confidently lead to the development of more effective immunotherapeutic approaches [97].
Here, we extensively analyzed the gene expression and genomic landscape of skin melanoma in the context of intratumoral immune cytolytic activity. We stratified primary and metastatic tumors according to a validated cytolytic gene expression signature, and discriminated a subgroup of them with noticeable T-cell reactivity. We show that CYT-high skin melanomas are significantly enriched for immune-related gene sets associated with activated CD8 + T cells, B cells, M1 macrophages, activated dendritic cells and NK cells, among others, corroborating the existence of an inflamed tumor microenvironment in these patients. In contrast, CYT-low skin melanomas were enriched for non-immune-related gene sets, primarily associated with monocytes and CD4 + T cells.
Abundant CD8 + T cell infiltrates are well known to exist in inflamed metastatic melanoma and drive the upregulation of PD-L1, IDO and Tregs in the tumor microenvironment [98–100]. Additionally, CD8 + T-cell infiltration associates with a better response of cancer patients to chemotherapy [101], neoadjuvant therapy [102] and anti–PD-1 immunotherapy [103]. A key role for tumor-infiltrating B cells in modulating the anti-tumor immune response, was recently proposed in skin melanoma [104]. B cells play an important prognostic role and can predict non-response to an immune checkpoint inhibitor in metastatic melanoma [104]. The enrichment of CYT-high tumors in M1 macrophages is in line with their suppressive role in cancer progression, while that of M2 macrophages in CYT-low tumors, underscores their role in favoring tumor growth and dissemination [105]. Activated dendritic cells are also important in the immune response against cancer cells and can be used as a strong independent prognostic factor [106]; whereas, activated NK cells can also efficiently kill malignant melanoma cells [107].
Taken together, these findings suggest that the stratification of skin melanoma patients according to their gene expression profiling can distinguish them between those having a strong cytolytic T-cell response, and those with a less privileged immune microenvironment. A similar, CYT-based stratification was previously performed for other types of cancer [14, 35, 36, 108].
We also showed that the mutational burden is higher among metastatic tumors, and that the increased mutational and neoantigen load that elicits an immune response by TILs, associates with high cytolytic levels in primary skin melanomas. In a similar study, a high mutational burden was previously associated with response in melanoma patients treated with immune checkpoint inhibition, but not with overall survival [109].
Both cytolytic subgroups of metastatic melanomas included driver mutations in some well-defined genes (BRAF, NRAS, CDKN2A, TP53 and PTEN) [80]. In addition, our data provide evidence that distinct cytolytic subgroups in primary and metastatic skin melanoma have different patterns of significantly mutated genes. We found that CYT-low primary tumors associate with mutations in BRAF, DMRT3, GSTA5, TP53, CRYBA4, STRA13, PRAMEF12, RNF32, SLC1A6, TEKT2 and NRAS; whereas, CYT-high primary tumors with LCK, GSTSF1L, BMF, CSTL1, DRGX, ENTPD3, CAMK4, GPR151 and LTF.
Similarly, we show that CYT-low metastatic tumors associate with mutations in GALNTL5, RGPD4, DNAJC5B, DEFA3, IRF2, COL2A1, VEGFC, WFDC11, PROL1, CD80, ADAM18, APOBEC3H and CSHL1; whereas, CYT-high metastatic tumors with PPP6C, DSG1, DEFB112, SPATA16, STARD6, RARRES2 PDE1A, IQCF3, KLK8, CDKN2A, SPANXN5, C16orf90, CD48, TXNDC3 and RAC1.
Although many of these mutations have been previously characterized in skin melanoma (e.g., BRAF, NRAS, TP53, PPP6C, CDKN2A, RAC1, etc.) [10, 80], this is the first report that appreciates them in association with tumors of different cytolytic activity.
We also provide evidence that different genomic structural variations are implicated in the progression of distinct cytolytic subgroups in skin melanoma. We show that CYT-low skin melanomas have a higher number of amplifications and deletions in their genome. These, include mainly recurrent NOTCH2 amplifications, and non-silent mutations and/or deletions in PD-L1, PD-L2, CDKN2A/B, PAX5, ETS1, BCL7, RAD51, JAK2, APAF1, FOXO3, CTNNA1 and IGF1, among others. These genomic structural variations have a profound effect on immune activation in these tumors.
Other distinct chromosomal aberrations included JAK2, KCNH2, MED21 and R3HDM2 amplifications or NRAS, RHOC, CDKN2B, ATM, CHEK1 and ETS1 deletions, and were associated with CYT-high tumors.
Many of them were previously appreciated in the genome of skin melanomas (e.g., MITF, NOTCH2, PD-L1 and JAK2 amplifications, or CDKN2A deletions [10]); however, the CNAs in genomic locations harboring PAX5, ETS1, BCL7, RAD51, APAF1, FOXO3 and CTNNA1 are novel.
Following these observations, we detected marked differences in the chromothriptic pattern between CYT-high and CYT-low tumors. These regarded the prevalence of chromothriptic events, the number of chromothriptic chromosomes per tumor and the cancer genes that they harbor, or the involvement of telomeres and centromeres in the chromothriptic regions. Such differences suggest that distinct mechanisms could give rise to chromothripsis across the two cytolytic subgroups in skin melanoma.
In addition, we highlight several targetable immune-related genes, other than PD-1, PD-L1/2 and CTLA-4, which may coordinately contribute to immunosuppression. IDO1 and LAG3 overexpression, which was evident in CYT-high skin melanomas, provides an alternative immunosuppressive barrier that is needed to hamper the antitumor activity of CTL and NK cells, in these tumors. This is something that was initially proposed in skin melanoma cell lines [110], and later corroborated at the human level, by us and others [14, 36]. We suggest that combined immune inhibition of these markers along with the PD-1/PD-L1 and/or the CTLA-4 axis, could provide a better therapeutic outcome for these patients [19].
Together with the synergistic effect of GZMA and PRF1 in patient survival, we provide evidence that CYT-high skin melanoma patients who received CTLA-4 or PD1/PD-L1/PD-L2 checkpoint blockade therapy (or a combination of both), have a markedly higher immunophenoscore and consequently, an expected clinical benefit, compared to CYT-low patients who did not receive checkpoint inhibition.
In the aggregate, we highlight new links between certain expressional or genomic changes and the activation of an antitumor immune response in skin melanoma. Our findings also corroborate that combination therapies with PD-1, PD-L1/2 and/or CTLA-4 blockade, should be able to overcome resistance and broaden the clinical benefit for these patients.
TCGA, The Cancer Genome Atlas
SKCM, skin melanoma
CYT, immune cytolytic activity
GZMA, granzyme A
PRF1, perforin-1
CTLA-4, cytotoxic T-lymphocyte-associated protein 4
PD-L1, Programmed death-ligand 1
PD-1, Programmed cell death protein 1
NK cells, natural killer cells
CTL, Cytotoxic T cells
SMGs, significantly mutated genes
VAF, variant allele frequency
MATH, mutant-allele tumor heterogeneity
MSigDB, Molecular Signatures Database
CNA, copy number alteration
CDN, classically defined neoepitope
ADN, alternatively defined neoepitope
Ethics approval and consent to participate: All data were retrieved from TCGA’s open access, therefore, no informed consent was necessary.
Consent for publication: Not applicable.
Availability of data and materials: All data generated or analyzed during this study are included in this published article [and its supplementary information files].
Competing interests: The authors declare that they have no competing interests.
Funding: No funding.
Author Contributions: AZ and CR developed methodology. CR, IGS and AZ acquired data, analyzed and interpreted them. AZ designed and supervised the study. AZ also wrote the paper.
Acknowledgements: We are grateful to Mr. Ioannis Apostolides for providing technical assistance with the antigen.garnish R package, as well as the contributors of the TCGA-SKCM dataset and the corresponding skin melanoma patients whose genetic and clinical data were analyzed.
- Bittner M, et al. Molecular classification of cutaneous malignant melanoma by gene expression profiling. Nature. 2000;406(6795):536–40.
- Brash DE, et al. A role for sunlight in skin cancer: UV-induced p53 mutations in squamous cell carcinoma. Proc Natl Acad Sci U S A. 1991;88(22):10124–8.
- Davies H, et al. Mutations of the BRAF gene in human cancer. Nature. 2002;417(6892):949–54.
- Huang FW, et al. Highly recurrent TERT promoter mutations in human melanoma. Science. 2013;339(6122):957–9.
- Leiter U, Garbe C. Epidemiology of melanoma and nonmelanoma skin cancer–the role of sunlight. Adv Exp Med Biol. 2008;624:89–103.
- Lin WM, et al. Modeling genomic diversity and tumor dependency in malignant melanoma. Cancer Res. 2008;68(3):664–73.
- Ruiz R, et al., The RhoJ-BAD signaling network: An Achilles' heel for BRAF mutant melanomas. 2017. 13(7): p. e1006913.
- Winsey SL, et al. A variant within the DNA repair gene XRCC3 is associated with the development of melanoma skin cancer. Cancer Res. 2000;60(20):5612–6.
- León-Letelier RA, Bonifaz LC, Fuentes-Pananá EM. OMIC signatures to understand cancer immunosurveillance and immunoediting: Melanoma and immune cells interplay in immunotherapy. J Leukoc Biol. 2019;105(5):915–33.
- Genomic Classification of Cutaneous Melanoma. Cell, 2015. 161(7): p. 1681–96.
- Danhier F, Feron O, Préat V. To exploit the tumor microenvironment: Passive and active tumor targeting of nanocarriers for anti-cancer drug delivery. J Control Release. 2010;148(2):135–46.
- Fu Q, et al. Prognostic value of tumor-infiltrating lymphocytes in melanoma: a systematic review and meta-analysis. Oncoimmunology. 2019;8(7):1593806.
- Kluger HM, et al. Characterization of PD-L1 Expression and Associated T-cell Infiltrates in Metastatic Melanoma Samples from Variable Anatomic Sites. Clin Cancer Res. 2015;21(13):3052–60.
- Rooney MS, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48–61.
- Mahoney KM, Freeman GJ, McDermott DF. The Next Immune-Checkpoint Inhibitors: PD-1/PD-L1 Blockade in Melanoma. Clin Ther. 2015;37(4):764–82.
- Hodi FS, et al. Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med. 2010;363(8):711–23.
- Luke JJ, et al. Targeted agents or immuno-oncology therapies as first-line therapy for BRAF-mutated metastatic melanoma: a real-world study. Future Oncol. 2019;15(25):2933–42.
- Topalian SL, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366(26):2443–54.
- Christodoulou MI, Zaravinos A. New Clinical Approaches and Emerging Evidence on Immune-Checkpoint Inhibitors as Anti-Cancer Therapeutics: CTLA-4 and PD-1 Pathways and Beyond. Crit Rev Immunol. 2019;39(5):379–408.
- Beatty GL, et al. First-in-Human Phase I Study of the Oral Inhibitor of Indoleamine 2,3-Dioxygenase-1 Epacadostat (INCB024360) in Patients with Advanced Solid Malignancies. Clin Cancer Res. 2017;23(13):3269–76.
- Rasku MA, et al. Transient T cell depletion causes regression of melanoma metastases. J Transl Med. 2008;6:12.
- Rech AJ, Vonderheide RH. Clinical use of anti-CD25 antibody daclizumab to enhance immune responses to tumor antigen vaccination by targeting regulatory T cells. Ann N Y Acad Sci. 2009;1174:99–106.
- Melero I, et al. Agonist antibodies to TNFR molecules that costimulate T and NK cells. Clin Cancer Res. 2013;19(5):1044–53.
- Moran AE, Kovacsovics-Bankowski M, Weinberg AD. The TNFRs OX40, 4-1BB, and CD40 as targets for cancer immunotherapy. Curr Opin Immunol. 2013;25(2):230–7.
- Nirschl CJ, Drake CG. Molecular pathways: coexpression of immune checkpoint molecules: signaling pathways and implications for cancer immunotherapy. Clin Cancer Res. 2013;19(18):4917–24.
- Nicholas C, Lesinski GB. Immunomodulatory cytokines as therapeutic agents for melanoma. Immunotherapy. 2011;3(5):673–90.
- Ngiow SF, et al. Anti-TIM3 antibody promotes T cell IFN-γ-mediated antitumor immunity and suppresses established tumors. Cancer Res. 2011;71(10):3540–51.
- John LB, et al. Anti-PD-1 antibody therapy potently enhances the eradication of established tumors by gene-modified T cells. Clin Cancer Res. 2013;19(20):5636–46.
- John LB, et al. Oncolytic virus and anti-4-1BB combination therapy elicits strong antitumor immunity against established cancer. Cancer Res. 2012;72(7):1651–60.
- Woo SR, et al. Immune inhibitory molecules LAG-3 and PD-1 synergistically regulate T-cell function to promote tumoral immune escape. Cancer Res. 2012;72(4):917–27.
- Takeda K, Okumura K, Smyth MJ. Combination antibody-based cancer immunotherapy. Cancer Sci. 2007;98(9):1297–302.
- Christofi T, et al., Current Perspectives in Cancer Immunotherapy. 2019. 11(10).
- Gide TN, et al. Primary and Acquired Resistance to Immune Checkpoint Inhibitors in Metastatic Melanoma. Clin Cancer Res. 2018;24(6):1260–70.
- Mermel CH, et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41.
- Zaravinos A, et al. Cytolytic activity correlates with the mutational burden and deregulated expression of immune checkpoints in colorectal cancer. J Exp Clin Cancer Res. 2019;38(1):364.
- Roufas C, et al. The Expression and Prognostic Impact of Immune Cytolytic Activity-Related Markers in Human Malignancies: A Comprehensive Meta-analysis. Front Oncol. 2018;8:27.
- Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
- Law CW, et al. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15(2):R29.
- Law CW, Alhamdoosh M. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. 2016. 5.
- Uhlén M, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.
- Thul PJ, Åkesson L. A subcellular map of the human proteome. 2017. 356(6340).
- Uhlen M, Zhang C. A pathology atlas of the human cancer transcriptome. 2017. 357(6352).
- Finotello F, et al., Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. 2019. 11(1): p. 34.
- Newman AM, Liu CL, Green MR. Robust enumeration of cell subsets from tissue expression profiles. 2015;12(5):453–7.
- Carter SL, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30(5):413–21.
- Lawrence MS, et al. Discovery and saturation analysis of cancer genes across 21 tumour types. Nature. 2014;505(7484):495–501.
- Lawrence MS, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214–8.
- Beroukhim R, et al. The landscape of somatic copy-number alteration across human cancers. Nature. 2010;463(7283):899–905.
- Yang J, et al. CTLPScanner: a web server for chromothripsis-like pattern detection. Nucleic Acids Res. 2016;44(W1):W252-8.
- Olshen AB, et al. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics. 2004;5(4):557–72.
- Forbes SA, et al. COSMIC: exploring the world's knowledge of somatic mutations in human cancer. Nucleic Acids Res. 2015;43(Database issue):D805-11.
- Campbell PJ. Telomeres and cancer: from crisis to stability to crisis to stability. Cell. 2012;148(4):633–5.
- Maciejowski J, de Lange T. Telomeres in cancer: tumour suppression and genome instability. Nat Rev Mol Cell Biol. 2017;18(3):175–86.
- Maciejowski J, et al. Chromothripsis and Kataegis Induced by Telomere Crisis. Cell. 2015;163(7):1641–54.
- Carreno BM, et al. Cancer immunotherapy. A dendritic cell vaccine increases the breadth and diversity of melanoma neoantigen-specific T cells. Science. 2015;348(6236):803–8.
- Gros A, et al. Prospective identification of neoantigen-specific lymphocytes in the peripheral blood of melanoma patients. Nat Med. 2016;22(4):433–8.
- Charoentong P, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017;18(1):248–62.
- Alexandrov LB, et al. Signatures of mutational processes in human cancer. Nature. 2013;500(7463):415–21.
- Díaz-Gay M, et al., Mutational Signatures in Cancer (MuSiCa): a web application to implement mutational signatures analysis in cancer samples. 2018. 19(1): p. 224.
- Nik-Zainal S, et al. Mutational processes molding the genomes of 21 breast cancers. Cell. 2012;149(5):979–93.
- Lada AG, et al., AID/APOBEC cytosine deaminase induces genome-wide kataegis. Biol Direct, 2012. 7: p. 47; discussion 47.
- Amelio I, et al., SynTarget: an online tool to test the synergetic effect of genes on survival outcome in cancer. 2016. 23(5): p. 912.
- Bolli N, et al. Heterogeneity of genomic evolution and mutational profiles in multiple myeloma. Nat Commun. 2014;5:2997.
- Chen X, et al. Recurrent somatic structural variations contribute to tumorigenesis in pediatric osteosarcoma. Cell Rep. 2014;7(1):104–12.
- Sakofsky CJ, et al. Break-induced replication is a source of mutation clusters underlying kataegis. Cell Rep. 2014;7(5):1640–8.
- Taylor BJ, et al. DNA deaminases induce break-associated mutation showers with implication of APOBEC3B and 3A in breast cancer kataegis. Elife. 2013;2:e00534.
- Boichard A, Tsigelny IF, Kurzrock R. High expression of PD-1 ligands is associated with kataegis mutational signature and APOBEC3 alterations. Oncoimmunology. 2017;6(3):e1284719.
- Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
- Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
- Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
- Ikeda H, Old LJ, Schreiber RD. The roles of IFN gamma in protection against tumor development and cancer immunoediting. Cytokine Growth Factor Rev. 2002;13(2):95–109.
- Ayers M, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017;127(8):2930–40.
- Han SB, et al. Pectic polysaccharide isolated from Angelica gigas Nakai inhibits melanoma cell metastasis and growth by directly preventing cell adhesion and activating host immune functions. Cancer Lett. 2006;243(2):264–73.
- Kepp O, et al. Extracellular nucleosides nucleotides as immunomodulators. 2017;280(1):83–92.
- Waugh KA, Leach SM, Slansky JE. Tolerance of Tumor-Specific T cells in Melanoma Metastases. J Clin Cell Immunol, 2016. 7(2).
- Bindea G, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39(4):782–95.
- Girardi M, et al. Regulation of cutaneous malignancy by gammadelta T cells. Science. 2001;294(5542):605–9.
- Ueno H. T follicular helper cells in human autoimmunity. Curr Opin Immunol. 2016;43:24–31.
- Fujimura T, Kambayashi Y, Aiba S. Crosstalk between regulatory T cells (Tregs) and myeloid derived suppressor cells (MDSCs) during melanoma growth. Oncoimmunology. 2012;1(8):1433–4.
- Hodis E, et al. A landscape of driver mutations in melanoma. Cell. 2012;150(2):251–63.
- Dabas N, et al. Diagnostic role of chromosomal instability in melanoma. J Skin Cancer. 2012;2012:914267.
- Hirsch D, et al. Chromothripsis and focal copy number alterations determine poor outcome in malignant melanoma. Cancer Res. 2013;73(5):1454–60.
- Shain AH, et al. The Genetic Evolution of Melanoma from Precursor Lesions. N Engl J Med. 2015;373(20):1926–36.
- Stephens PJ, et al. Massive genomic rearrangement acquired in a single catastrophic event during cancer development. Cell. 2011;144(1):27–40.
- Brás A, Rodrigues AS, Rueff J. Copy number variations and constitutional chromothripsis (Review). Biomed Rep. 2020;13(3):11.
- Maciejowski J, et al., APOBEC3-dependent kataegis and TREX1-driven chromothripsis during telomere crisis. 2020.
- Ernst A, et al. Telomere dysfunction and chromothripsis. Int J Cancer. 2016;138(12):2905–14.
- Zhang CZ, et al. Chromothripsis from DNA damage in micronuclei. Nature. 2015;522(7555):179–84.
- Ly P, et al., Selective Y centromere inactivation triggers chromosome shattering in micronuclei and repair by non-homologous end joining. 2017. 19(1): p. 68–75.
- Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science. 2015;348(6230):69–74.
- McGranahan N, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science. 2016;351(6280):1463–9.
- Yarchoan M, et al. Targeting neoantigens to augment antitumour immunity. Nat Rev Cancer. 2017;17(4):209–22.
- Rizvi NA, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124–8.
- Spencer KR, et al. Biomarkers for Immunotherapy: Current Developments and Challenges. Am Soc Clin Oncol Educ Book. 2016;35:e493–503.
- Van Allen EM, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350(6257):207–11.
- Hugo W, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2017;168(3):542.
- Dunn GP, Old LJ, Schreiber RD. The immunobiology of cancer immunosurveillance and immunoediting. Immunity. 2004;21(2):137–48.
- Gajewski TF, Louahed J, Brichard VG. Gene signature in melanoma associated with clinical activity: a potential clue to unlock cancer immunotherapy. Cancer J. 2010;16(4):399–403.
- Harlin H, et al. Chemokine expression in melanoma metastases associated with CD8 + T-cell recruitment. Cancer Res. 2009;69(7):3077–85.
- Spranger S, et al. Up-regulation of PD-L1, IDO, and T(regs) in the melanoma tumor microenvironment is driven by CD8(+) T cells. Sci Transl Med. 2013;5(200):200ra116.
- Fridman WH, et al. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12(4):298–306.
- Denkert C, et al. Tumor-associated lymphocytes as an independent predictor of response to neoadjuvant chemotherapy in breast cancer. J Clin Oncol. 2010;28(1):105–13.
- Tumeh PC, et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature. 2014;515(7528):568–71.
- Selitsky SR, et al. Prognostic value of B cells in cutaneous melanoma. Genome Med. 2019;11(1):36.
- Falleni M, et al. M1 and M2 macrophages' clinicopathological significance in cutaneous melanoma. Melanoma Res. 2017;27(3):200–10.
- Ladányi A, et al. Density of DC-LAMP(+) mature dendritic cells in combination with activated T lymphocytes infiltrating primary cutaneous melanoma is a strong independent prognostic factor. Cancer Immunol Immunother. 2007;56(9):1459–69.
- Pietra G, et al. Natural killer cells kill human melanoma cells with characteristics of cancer stem cells. Int Immunol. 2009;21(7):793–801.
- Balli D, et al. Immune Cytolytic Activity Stratifies Molecular Subsets of Human Pancreatic Cancer. Clin Cancer Res. 2017;23(12):3129–38.
- Morrison C, et al. Predicting response to checkpoint inhibitors in melanoma beyond PD-L1 and mutational burden. J Immunother Cancer. 2018;6(1):32.
- Pietra G, et al. Melanoma cells inhibit natural killer cell function by modulating the expression of activating receptors and cytolytic activity. Cancer Res. 2012;72(6):1407–15.
This is a list of supplementary files associated with this preprint. Click to download.
- TableS4.pdf
Table S4. Cancer genes within chromosomal regions that harbored chromothripsis, with events of genomic gains or losses, across cytolytic high (CYT-high) or low (CYT-low) skin melanomas.
- TableS4.pdf
Table S4. Cancer genes within chromosomal regions that harbored chromothripsis, with events of genomic gains or losses, across cytolytic high (CYT-high) or low (CYT-low) skin melanomas.
- TableS4.pdf
Table S4. Cancer genes within chromosomal regions that harbored chromothripsis, with events of genomic gains or losses, across cytolytic high (CYT-high) or low (CYT-low) skin melanomas.
- TableS3.pdf
Table S3. Summary table of the key statistics for the top differentially expressed genes in primary and metastatic skin melanoma, compared to the normal skin. logFC: estimate of the log2-fold-change corresponding to the effect or contrast; AveExpr: average log2-expression; t: moderated t-statistic; F: moderated F-statistic; P.Value: raw p-value; adj.p: adjusted p-value; B: log-odds that the gene is differentially expressed.
- TableS3.pdf
Table S3. Summary table of the key statistics for the top differentially expressed genes in primary and metastatic skin melanoma, compared to the normal skin. logFC: estimate of the log2-fold-change corresponding to the effect or contrast; AveExpr: average log2-expression; t: moderated t-statistic; F: moderated F-statistic; P.Value: raw p-value; adj.p: adjusted p-value; B: log-odds that the gene is differentially expressed.
- TableS3.pdf
Table S3. Summary table of the key statistics for the top differentially expressed genes in primary and metastatic skin melanoma, compared to the normal skin. logFC: estimate of the log2-fold-change corresponding to the effect or contrast; AveExpr: average log2-expression; t: moderated t-statistic; F: moderated F-statistic; P.Value: raw p-value; adj.p: adjusted p-value; B: log-odds that the gene is differentially expressed.
- Suppl.Table2.xlsx
Table S2. Antibody staining, intensity, quantity and location of GZMA and PRF1 proteins across 12 skin melanomas. Immunohistochemistry data were extracted from the Human Protein Atlas using antibody-based protein profiling. The numbers indicate the skin melanoma patients expressing each gene.
- Suppl.Table2.xlsx
Table S2. Antibody staining, intensity, quantity and location of GZMA and PRF1 proteins across 12 skin melanomas. Immunohistochemistry data were extracted from the Human Protein Atlas using antibody-based protein profiling. The numbers indicate the skin melanoma patients expressing each gene.
- Suppl.Table2.xlsx
Table S2. Antibody staining, intensity, quantity and location of GZMA and PRF1 proteins across 12 skin melanomas. Immunohistochemistry data were extracted from the Human Protein Atlas using antibody-based protein profiling. The numbers indicate the skin melanoma patients expressing each gene.
- TableS1.pdf
Table S1. Clinical information of the patients in the TCGA-SKCM dataset.
- TableS1.pdf
Table S1. Clinical information of the patients in the TCGA-SKCM dataset.
- TableS1.pdf
Table S1. Clinical information of the patients in the TCGA-SKCM dataset.
- SupplFigure8.pdf
Figure S8. BRAF, NRAS, TP53 and NF1 mutations across the SKCM dataset and association with CYT, showing no preference for these four gene mutations in any cytolytic subgroup.
- SupplFigure8.pdf
Figure S8. BRAF, NRAS, TP53 and NF1 mutations across the SKCM dataset and association with CYT, showing no preference for these four gene mutations in any cytolytic subgroup.
- SupplFigure8.pdf
Figure S8. BRAF, NRAS, TP53 and NF1 mutations across the SKCM dataset and association with CYT, showing no preference for these four gene mutations in any cytolytic subgroup.
- SupplFigure7.pdf
Figure S7. 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.
- SupplFigure7.pdf
Figure S7. 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.
- SupplFigure7.pdf
Figure S7. 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.
- Suppl.Figure6.pdf
Figure S6. The expression of the immunoinhibitor LAG3, immunostimulators (CD70, LTA, NT5E and ENTPD1), markers for activated CD8+ T cells (NKG7, CD3E, GZMA, GZMH, GZMK), MDSC (CD2), activated dendritic cells (UBD, C1QB, C1QC), NK cells (FASLG and FAS) and chemokines the CXCL9, CXCL10, CXCL11 and CXCL13, was significantly higher across CYT-high (primary and metastatic) cutaneous melanomas.
- Suppl.Figure6.pdf
Figure S6. The expression of the immunoinhibitor LAG3, immunostimulators (CD70, LTA, NT5E and ENTPD1), markers for activated CD8+ T cells (NKG7, CD3E, GZMA, GZMH, GZMK), MDSC (CD2), activated dendritic cells (UBD, C1QB, C1QC), NK cells (FASLG and FAS) and chemokines the CXCL9, CXCL10, CXCL11 and CXCL13, was significantly higher across CYT-high (primary and metastatic) cutaneous melanomas.
- Suppl.Figure6.pdf
Figure S6. The expression of the immunoinhibitor LAG3, immunostimulators (CD70, LTA, NT5E and ENTPD1), markers for activated CD8+ T cells (NKG7, CD3E, GZMA, GZMH, GZMK), MDSC (CD2), activated dendritic cells (UBD, C1QB, C1QC), NK cells (FASLG and FAS) and chemokines the CXCL9, CXCL10, CXCL11 and CXCL13, was significantly higher across CYT-high (primary and metastatic) cutaneous melanomas.
- Suppl.Figure5.pdf
Figure S5. The unsupervised hierarchical clustering of the two cytolytic subgroups using the two IFN-γ-related gene signatures proposed by Ayers, et al. (2017). These are the “IFN-γ signature” composed of IDO1, CXCL10, CXCL9, HLA-DRA, STAT1 and IFNG (a), and the “expanded IFN-γ signature” composed of CD30, IDO1, CIITA, CD3E, CCL5, GZMK, CD2, HLA-DRA, CXCL13, IL2RG, NKG7, HLA-E, CXCR6, LAG3, TAGAP, CXCL10, STAT1 and GZMB (b). Both signatures shows are enriched in CYT-high skin melanomas. (c) Boxplots depict the individual gene expression levels of each signature gene in the two cytolytic subsets of skin melanoma (SKCM).
- Suppl.Figure5.pdf
Figure S5. The unsupervised hierarchical clustering of the two cytolytic subgroups using the two IFN-γ-related gene signatures proposed by Ayers, et al. (2017). These are the “IFN-γ signature” composed of IDO1, CXCL10, CXCL9, HLA-DRA, STAT1 and IFNG (a), and the “expanded IFN-γ signature” composed of CD30, IDO1, CIITA, CD3E, CCL5, GZMK, CD2, HLA-DRA, CXCL13, IL2RG, NKG7, HLA-E, CXCR6, LAG3, TAGAP, CXCL10, STAT1 and GZMB (b). Both signatures shows are enriched in CYT-high skin melanomas. (c) Boxplots depict the individual gene expression levels of each signature gene in the two cytolytic subsets of skin melanoma (SKCM).
- Suppl.Figure5.pdf
Figure S5. The unsupervised hierarchical clustering of the two cytolytic subgroups using the two IFN-γ-related gene signatures proposed by Ayers, et al. (2017). These are the “IFN-γ signature” composed of IDO1, CXCL10, CXCL9, HLA-DRA, STAT1 and IFNG (a), and the “expanded IFN-γ signature” composed of CD30, IDO1, CIITA, CD3E, CCL5, GZMK, CD2, HLA-DRA, CXCL13, IL2RG, NKG7, HLA-E, CXCR6, LAG3, TAGAP, CXCL10, STAT1 and GZMB (b). Both signatures shows are enriched in CYT-high skin melanomas. (c) Boxplots depict the individual gene expression levels of each signature gene in the two cytolytic subsets of skin melanoma (SKCM).
- SupplFigure4.pdf
Figure S4. Purity (a) and ploidy (b) scores across the two cytolytic subsets of primary and metastatic skin melanomas, using ABSOLUTE.
- SupplFigure4.pdf
Figure S4. Purity (a) and ploidy (b) scores across the two cytolytic subsets of primary and metastatic skin melanomas, using ABSOLUTE.
- SupplFigure4.pdf
Figure S4. Purity (a) and ploidy (b) scores across the two cytolytic subsets of primary and metastatic skin melanomas, using ABSOLUTE.
- SupplFigure3.pdf
Figure S3. (a) The plots depict the voom, mean-variance trend for primary and metastatic skin melanomas, respectively. (b) Principal-components analysis (PCA) for the two cytolytic subgroups of primary and metastatic skin melanomas. (c) Volcano plot depicting the significantly differentially expressed genes in primary (left) and metastatic (right) skin melanomas. The top 50 differentially upregulated genes are highlighted in blue. (d) The Venn-diagram depicts the unique differential genes within primary (13 genes) and metastatic skin melanomas (8,121 genes), as well as the co-deregulated genes (1,054 genes) between the two tumor types.
- SupplFigure3.pdf
Figure S3. (a) The plots depict the voom, mean-variance trend for primary and metastatic skin melanomas, respectively. (b) Principal-components analysis (PCA) for the two cytolytic subgroups of primary and metastatic skin melanomas. (c) Volcano plot depicting the significantly differentially expressed genes in primary (left) and metastatic (right) skin melanomas. The top 50 differentially upregulated genes are highlighted in blue. (d) The Venn-diagram depicts the unique differential genes within primary (13 genes) and metastatic skin melanomas (8,121 genes), as well as the co-deregulated genes (1,054 genes) between the two tumor types.
- SupplFigure3.pdf
Figure S3. (a) The plots depict the voom, mean-variance trend for primary and metastatic skin melanomas, respectively. (b) Principal-components analysis (PCA) for the two cytolytic subgroups of primary and metastatic skin melanomas. (c) Volcano plot depicting the significantly differentially expressed genes in primary (left) and metastatic (right) skin melanomas. The top 50 differentially upregulated genes are highlighted in blue. (d) The Venn-diagram depicts the unique differential genes within primary (13 genes) and metastatic skin melanomas (8,121 genes), as well as the co-deregulated genes (1,054 genes) between the two tumor types.
- SupplFigure2.pdf
Figure S2. Data normalization and filtering on the RNA-seq read counts. Read counts were converted to log2-counts-per-million (logCPM) and voom transformed using the limma package. (a) Un-normalized logCPM, (b) Voom transformed logCPM.
- SupplFigure2.pdf
Figure S2. Data normalization and filtering on the RNA-seq read counts. Read counts were converted to log2-counts-per-million (logCPM) and voom transformed using the limma package. (a) Un-normalized logCPM, (b) Voom transformed logCPM.
- SupplFigure2.pdf
Figure S2. Data normalization and filtering on the RNA-seq read counts. Read counts were converted to log2-counts-per-million (logCPM) and voom transformed using the limma package. (a) Un-normalized logCPM, (b) Voom transformed logCPM.
- Suppl.Figure1.pdf
Figure S1. (a, b) The reconstructed mutational profiles and COSMIC signatures contributions did not differ between CYT-high and CYT-low (primary and metastatic) SKCM tumors (C>T relative contribution). The residual sum of squares (RSS) is a measure of the efficiency of the original mutational profile reconstruction. The cosine similarity between profiles presents a direct measure of the correspondence between the two depicted profiles in a 0–1 range. Analysis was performed using MuSiCa. (c) Rainfall plots depicting the log10 of the intermutation distance versus genomic position for two representative skin melanoma samples with evidence of kataegis: one CYT-high (TCGA-3N-A9WC-06A-11D-A38G-08) and onw CYT-low (TCGA-D3-A8GL-06A-11D-A372-08) skin tumor. The dots represent the mutations and the different colors represent each of the 6 subtypes of base substitutions: C>A, C>G, C>T, T>A, T>C, and T>G. Arrowheads depict clusters of mutations in kataegis.
- Suppl.Figure1.pdf
Figure S1. (a, b) The reconstructed mutational profiles and COSMIC signatures contributions did not differ between CYT-high and CYT-low (primary and metastatic) SKCM tumors (C>T relative contribution). The residual sum of squares (RSS) is a measure of the efficiency of the original mutational profile reconstruction. The cosine similarity between profiles presents a direct measure of the correspondence between the two depicted profiles in a 0–1 range. Analysis was performed using MuSiCa. (c) Rainfall plots depicting the log10 of the intermutation distance versus genomic position for two representative skin melanoma samples with evidence of kataegis: one CYT-high (TCGA-3N-A9WC-06A-11D-A38G-08) and onw CYT-low (TCGA-D3-A8GL-06A-11D-A372-08) skin tumor. The dots represent the mutations and the different colors represent each of the 6 subtypes of base substitutions: C>A, C>G, C>T, T>A, T>C, and T>G. Arrowheads depict clusters of mutations in kataegis.
- Suppl.Figure1.pdf
Figure S1. (a, b) The reconstructed mutational profiles and COSMIC signatures contributions did not differ between CYT-high and CYT-low (primary and metastatic) SKCM tumors (C>T relative contribution). The residual sum of squares (RSS) is a measure of the efficiency of the original mutational profile reconstruction. The cosine similarity between profiles presents a direct measure of the correspondence between the two depicted profiles in a 0–1 range. Analysis was performed using MuSiCa. (c) Rainfall plots depicting the log10 of the intermutation distance versus genomic position for two representative skin melanoma samples with evidence of kataegis: one CYT-high (TCGA-3N-A9WC-06A-11D-A38G-08) and onw CYT-low (TCGA-D3-A8GL-06A-11D-A372-08) skin tumor. The dots represent the mutations and the different colors represent each of the 6 subtypes of base substitutions: C>A, C>G, C>T, T>A, T>C, and T>G. Arrowheads depict clusters of mutations in kataegis.
Badges
Prescreen
This manuscript passed our Prescreen assessment
Complete author information
Appropriate declaration statements
Potential risks to human health
Prescreen
History
CURRENT STATUS: Posted
Version 1
Posted 15 Sep, 2020
No community comments so far
Metrics
Comments: 0
PDF Downloads: ...
HTML Views: ...
Subject Areas
Cancer Biology
Learn more about our company.
This is a preprint. It has not completed peer review.
Research
Distinct Genomic Features Across Cytolytic Subgroups in Skin Melanoma
Apostolos Zaravinos
Department of Basic Medical Sciences, College of Medicine, Member of QU Health, Qatar University, 2713 Doha, Qatar.
Corresponding Author
ORCiD: https://orcid.org/0000-0003-4625-5562
Badges
Prescreen
This manuscript passed our Prescreen assessment
Complete author information
Appropriate declaration statements
Potential risks to human health
Prescreen
History
CURRENT STATUS: Posted
Version 1
Posted 15 Sep, 2020
No community comments so far
Metrics
Comments: 0
PDF Downloads: ...
HTML Views: ...
Subject Areas
Cancer Biology
License:
This work is licensed under a Creative Commons Attribution 4.0 International License. Read Full License
Background: Skin melanoma is a highly immunogenic cancer with extensive genetic and transcriptional diversity. Despite the enthusiastic clinical results seen in advanced-stage metastatic melanoma patients treated with immune checkpoint inhibition, a subgroup of them later relapse and develop acquired resistance. The intratumoral immune cytolytic activity (CYT) reflects the ability of cytotoxic T cells and NK cells to eliminate cancer cells, and is associated with improved patient survival. Here, we questioned whether CYT associates with different genomic profiles in skin melanoma.
Methods: We explored the TCGA-SKCM dataset and stratified patients to distinct subgroups of cytolytic activity. The tumor immune contexture, somatic mutations and recurrent copy number aberrations were calculated using quanTIseq, MutSigCV and GISTIC2. Chromothriptic events were explored using CTLPScanner and cancer neoepitopes were predicted with antigen.garnish. Each tumor's immunophenoscore was calculated using Immunophenogram.
Results: Metastatic skin melanomas had significantly higher CYT compared to primary tumors. We also assessed enrichment for immune-related gene sets within CYT-high tumors; whereas, CYT-low tumors were enriched for non-immune related genes sets. In addition, distinct mutational and neoantigenic loads, primarily composed of C>T transitions, along with specific types of copy number aberrations, characterized each cytolytic subgroup. We found a broader pattern of chromothripsis across CYT-low tumors, where non-telomeric chromosomal regions harboring chromothripsis, contained a higher number of cancer genes. CYT-high patients had markedly higher immunophenoscore and should consequently, display an expected clinical benefit compared to CYT-low patients who either received or not, checkpoint inhibition.
Conclusions: Overall, our data highlight the existence of distinct genomic features across cytolytic subgroups in skin melanoma.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Cutaneous melanoma is a very aggressive and highly immunogenic cancer with extensive genetic and transcriptional diversity [1–8]. The tumor usually results after exposure to UV radiation, and has the highest rate of somatic mutations and neoantigens among all cancer types [9]. The Cancer Genome Atlas (TGCA) Network recently classified skin melanoma into mutant BRAF, mutant RAS, mutant NF1, and Triple-WT genetic subtypes, based on the pattern of the top mutated genes [10].
The tumor’s microenvironment comprises a very heterogeneous cell population, including fibroblasts, lymphocytes, macrophages and other immune cells, as well as, adipocytes and cells that form the structural elements of skin blood vessels dashed in the extracellular matrix [11]. A high content of tumor-infiltrating lymphocytes (TILs) in the tumor’s microenvironment associates with a favorable prognosis and improved overall survival of the patients [12, 13]. An increased immune cytolytic activity (CYT), defined by the attempt of the cytotoxic T cells and natural killer (NK) cells to eliminate cancer cells via the secretion of GZMA and PRF1, has also been associated with improved patient survival [14].
With the advent of immunotherapies, the tumor’s management has shifted from cytokine-based treatment to immune checkpoint inhibition, primarily of cytotoxic T-lymphocyte-associated antigen-4 (CTLA-4) and programmed cell-death protein 1 (PD-1) or its ligands (PD-L1 and PD-L2) [15–18]. Monoclonal antibodies including Ipilimumab, Tremelimumab (anti-CTLA-4), Pembrolizumab, Nivolumab, Cemiplimab (anti-PD-1), Atezolizumab, Avelumab and Durvalumab (anti-PD-L1), among others currently being on clinical trial [19], have enthusiastically provided long-lasting responses and improved survival in patients with advanced-stage metastatic melanoma [16].
Early clinical results have also been stated with IDO inhibitors [20] and Treg depletion targeting surface CD25 [21, 22]. Many other immunotherapeutic agents are currently under study, targeting a range of positive or negative immunoregulatory molecules within the tumor microenvironment. These, include agonistic antibodies against 4-1BB, Ox40, Inducible T-cell COStimulator (ICOS), and CD40; blocking antibodies against Lymphocyte-activation gene 3 (LAG3), B7-H3, B7-H4, Tim3, and killer inhibitory receptors (KIRs); as well as cytokines IL-7, IL-15 and IL-21 [19, 23–32], among others.
Due to the complexity of immune regulation in vivo, combination immunotherapies are expected to provide a better therapeutic benefit for the patients [19]. Nevertheless, a subgroup of responders to immunotherapy, later relapse and develop acquired resistance; whereas, others do not respond at all (primary resistance) [33]. Therefore, a better understanding of the underlying resistance mechanisms is emergently needed.
Here, we questioned whether CYT associates with different genomic profiles in skin melanoma. We investigated how the immune landscape in these tumors relates to somatic mutations, mutational signatures, copy number aberrations, chromothriptic events, the expression of immune checkpoints or other immune-related markers, and the presence of different types of immune cells within the tumor microenvironment of primary or metastatic tumors. We also predicted how skin melanoma patients of different cytolytic activity respond to immune checkpoint blockade therapy. Overall, we provide enough evidence of the existence of distinct genomic features across cytolytic subgroups of skin melanoma.
Skin melanoma (SKCM) data extraction
We extracted the clinical information, Mutation Annotation files (MAF) and mRNA-Seq ‘level 3’ data of a total of 470 primary and metastatic skin melanomas (103 primary tumors, 68 distant metastasis, 74 regional cutaneous or subcutaneous tissues and 222 regional lymph nodes), along with their matched peripheral blood, from the Cancer Genome Atlas’ (TCGA-SKCM dataset) GDC Data Portal (https://portal.gdc.cancer.gov/). Three samples were non-informative and were thus, excluded from subsequent analyses. In addition, we accessed each patient’s gene-level, zero-centered, focal copy-number data analyzed by GISTIC (v2.0.22) [34] from Broad GDAC Firehose (https://gdac.broadinstitute.org/).
Cytolytic Activity Calculation And Downstream Rna-seq Analysis
We calculated each SKCM patient’s levels of activity (CYT), as the geometric mean of GZMA and PRF1 [14, 35, 36]. GZMA leads to caspase-independent apoptosis, while PRF1 forms pores in the target tumor cells facilitating the entry of granzymes into them. Gene expression values were presented in Transcripts Per Million (TPM).
We then divided patients into two cohorts, the upper 25th quartile of the cytolytic index (CYT-high) and the lower 25th quartile (CYT-low), each with an identical admixture of histology-stage combinations. All subsequent comparisons were made between CYT-high and CYT-low (metastatic or primary) skin melanomas. P-values were adjusted using the Benjamini-Hochberg (BH) method.
We used gene set variation analysis (GSVA) to refine alterations in pathway activity [37] and clustered tumors hierarchically, with complete linkage as the distance metric. We compared gene expression between the two cytolytic subsets of melanomas, first by estimating the mean-variance relationship of the log-counts using voom [38] and then analyzing them with limma [39]. We considered differentially expressed genes, those having a BH-adjusted p-value < 0.1. Graphs were plotted with ggplot2.
Immunohistochemistry (IHC) data of GZMA and PRF1 protein expression were extracted from the Human Protein Atlas [40–42], and further analyzed. The antibodies used in IHC were as follows: rabbit pAb anti-GZMA, HPA054134, 1:200 dilution, Sigma-Aldrich, Atlas Antibodies Cat#HPA054134, RRID:AB_2682395; Antigen retrieval was performed using HIER pH6; and mouse mAb anti-PRF1, CAB002436, 1:10 dilution, Leica Biosystems, Cat#NCL-PERFORIN, RRID:AB_563955; Antigen retrieval was performed using HIER pH6.
Cell Type Fractions And Heterogeneity
We quantified the tumor immune contexture from the RNA-seq data, using quanTIseq (https://icbi.i-med.ac.at/software/quantiseq/doc/) [43]. Additionally, we used CIBERSORT to identify fractions of immune subpopulations in each cytolytic subgroup’s tissue (https://cibersort.stanford.edu/) [44]. Tumor purity and ploidy estimates were produced using ABSOLUTE [45] and plotted for each cytolytic subgroup of primary and metastatic tumors.
Detection Of Somatic Mutations And Copy Number Aberrations (scna)
We calculated the most significantly mutated genes (SMG) in each CYT subtype of primary and metastatic skin melanomas, using MutSigCV (v1.3.01) [46, 47] with an FDR (q-value) = 0.1 as threshold of significance. We further employed GISTIC (v2.0.22) [34, 48] to identify genomic regions with significant amplifications or deletions in each sample. A G-score considering the amplitude of the aberration, as well as the frequency of its occurrence across each tumor sample was assigned. We calculated the somatic copy number alterations in each sample by taking the sum of segment mean changes ≥ 0.6 and ≤ -0.4 between somatic and normal samples. Regions with FDR q-values < 0.1 were considered significant. Both mutational and copy number analyses were performed using the Broad Institute’s GenePattern platform (http://software.broadinstitute.org/cancer/software/genepattern).
Detection Of Chromothriptic Events
Chromothriptic events in CYT-high or low skin melanomas were investigated using CTLPScanner [49]. In brief, after downloading segmentation data (level 3) of single nucleotide polymorphism (SNP) arrays from the TCGA-SKCM dataset, we implemented DNA copy number segmentation data, using the circular binary segmentation algorithm [50], and the NCBI38/hg38 genome assembly. For the detection of chomothripsis or chomothripsis-like regions, we used copy number status change ≥ 20 times, log10 of likelihood ratio ≥ 8, minimum segment size = 10 kb, and signal distance between adjacent segments = 0.3. Shattered chromosomomal regions were visualized based on the signal value for genomic gains (≥ 0.15) or losses (≤ -0.15) and further highlighted. Finally, chromothripsis-located genes were annotated using the COSMIC database [51].
Consensus centromeric regions were extracted from http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/centromeres.txt.gz. Telomeric regions were defined as 20 kb from each side of the chromosomal ends. Intersection between chromothripsis events and centromeric or telomeric regions was defined as any overlap between the event coordinates and the region coordinates.
Histograms of copy number segment switches were drawn to identify events of chromothripsis that arise in conjunction with telomere crisis, i.e., a period of genome instability during tumorigenesis when depletion of the telomere reserve generates unstable dicentric chromosomes [52–54]. The patterns and distribution of copy number state switches indicate the emergence and temporal order of the major genomic rearrangement events.
Detection Of Cancer Neoepitopes And Immunophenoscores
We predicted the cancer neoepitopes using “antigen.garnish”, with the following MHC molecules, as previously described [35, 55, 56]: H-2-IAb, H-2-IAd, HLA-DPA1*01:03, HLA-DPA1*02:01, HLA-DPA1*02:01, HLA-DPA1*03:01, HLA-DPB1*03:01, HLA-DQA1*01:01, HLA-DQA1*01:02, HLA-DQA1*03:01, HLA-DQA1*04:01, HLA-DQA1*05:01, HLA-DQA1*05:01, HLA-DRB1*01:01. We classified the peptides that were predicted to bind MHC with high affinity (IC50 < 50 nM) or with greatly improved affinity compared to their wild-type counterparts, either as classically (CDNs) or alternatively defined neoepitopes (ADNs).
The Cancer Immunome Database (TCIA) (https://tcia.at/) was also queried for cellular composition of neoantigens among the two cytolytic subsets of skin melanomas [57]. The immunophenoscore (IPS) of each cytolytic subgroup’s patients was calculated using Immunophenogram, as previously described [57], and immunophenograms were constructed to visualize the immunophenotypes of each tumor sample. The IPS ranged from 0 to 10, according to the expression of MHC molecules, immunomodulators, effector cells (activated CD8 + T cells and CD4 + T cells, effector memory (Tem) CD8 + and CD4 + cells), and suppressor cells, including regulatory T cells (Tregs) and myeloid-derived suppressor cells (MDSCs) [57].
Mutational Signatures And Rainfall Plots
We extracted mutational signatures as previously described [58, 59] and compared them against the validated mutational signatures provided by COSMIC (http://cancer.sanger.ac.uk/cosmic/signatures). We further applied cosine similarity against COSMIC's validated signatures to identify the best matches within signatures. We reconstructed the mutational spectrum of CYT-high and CYT-low skin melanomas, using the residual sum of squares (RSS) to measure the efficiency of the reconstruction of the original mutational profile.
The genomic regions of each cytolytic subgroup that exhibited rainfall plots, i.e., contained ≥ 6 consecutive mutations with an average inter-mutation distance of ≤ 1 bp, were defined as kataegis [58, 60, 61].
Tumor Heterogeneity
Tumor heterogeneity in SKCM tumors was inferred by clustering variant allele frequencies (VAF). The extent of each skin melanoma’s intra-tumor heterogeneity was quantitatively calculated by the width of the VAF distribution. We assigned a mutant-allele tumor heterogeneity (MATH) score to each SKCM sample, as previously described [35].
Patient Survival And Synergistic Target Analysis
The overall survival of the melanoma patients was performed using data extracted from the human skin cutaneous melanoma dataset (TCGA-SKCM). Kaplan–Meier curves using log rank (Mantel Cox) test estimated the patient survival. We examined whether GZMA and PRF1 act synergistically on patient survival, using SynTarget [62].
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 [14]. 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 [67], 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 [68] producing a matrix of counts. The read counts were converted to log2-counts-per-million (logCPM) and voom transformed using limma [69] 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 [45]. 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) [70], 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 [71]. 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 [72]. 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 [72].
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) [73], ecto-5′-nucleotidase (NT5E) [74] and ectonucleoside triphosphate diphosphohydrolase 1 (ENTPD1) [75], 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) [76]. 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 [77], autoimmunity [78] and immunosuppression [79].
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 [80] 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 [45] 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 [62].
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 [84] and involves impaired DNA repair [85]. 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 [92].
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) [57], 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 [93], or clonal neoantigens [91], have all been proposed as markers; however, none of them has been fully validated, yet [94].
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 [95] and anti-PD-1 inhibitory molecules [96], 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 [95].
A better understanding of the immunobiology of cancer immunosurveillance and immunoediting in skin melanoma, will confidently lead to the development of more effective immunotherapeutic approaches [97].
Here, we extensively analyzed the gene expression and genomic landscape of skin melanoma in the context of intratumoral immune cytolytic activity. We stratified primary and metastatic tumors according to a validated cytolytic gene expression signature, and discriminated a subgroup of them with noticeable T-cell reactivity. We show that CYT-high skin melanomas are significantly enriched for immune-related gene sets associated with activated CD8 + T cells, B cells, M1 macrophages, activated dendritic cells and NK cells, among others, corroborating the existence of an inflamed tumor microenvironment in these patients. In contrast, CYT-low skin melanomas were enriched for non-immune-related gene sets, primarily associated with monocytes and CD4 + T cells.
Abundant CD8 + T cell infiltrates are well known to exist in inflamed metastatic melanoma and drive the upregulation of PD-L1, IDO and Tregs in the tumor microenvironment [98–100]. Additionally, CD8 + T-cell infiltration associates with a better response of cancer patients to chemotherapy [101], neoadjuvant therapy [102] and anti–PD-1 immunotherapy [103]. A key role for tumor-infiltrating B cells in modulating the anti-tumor immune response, was recently proposed in skin melanoma [104]. B cells play an important prognostic role and can predict non-response to an immune checkpoint inhibitor in metastatic melanoma [104]. The enrichment of CYT-high tumors in M1 macrophages is in line with their suppressive role in cancer progression, while that of M2 macrophages in CYT-low tumors, underscores their role in favoring tumor growth and dissemination [105]. Activated dendritic cells are also important in the immune response against cancer cells and can be used as a strong independent prognostic factor [106]; whereas, activated NK cells can also efficiently kill malignant melanoma cells [107].
Taken together, these findings suggest that the stratification of skin melanoma patients according to their gene expression profiling can distinguish them between those having a strong cytolytic T-cell response, and those with a less privileged immune microenvironment. A similar, CYT-based stratification was previously performed for other types of cancer [14, 35, 36, 108].
We also showed that the mutational burden is higher among metastatic tumors, and that the increased mutational and neoantigen load that elicits an immune response by TILs, associates with high cytolytic levels in primary skin melanomas. In a similar study, a high mutational burden was previously associated with response in melanoma patients treated with immune checkpoint inhibition, but not with overall survival [109].
Both cytolytic subgroups of metastatic melanomas included driver mutations in some well-defined genes (BRAF, NRAS, CDKN2A, TP53 and PTEN) [80]. In addition, our data provide evidence that distinct cytolytic subgroups in primary and metastatic skin melanoma have different patterns of significantly mutated genes. We found that CYT-low primary tumors associate with mutations in BRAF, DMRT3, GSTA5, TP53, CRYBA4, STRA13, PRAMEF12, RNF32, SLC1A6, TEKT2 and NRAS; whereas, CYT-high primary tumors with LCK, GSTSF1L, BMF, CSTL1, DRGX, ENTPD3, CAMK4, GPR151 and LTF.
Similarly, we show that CYT-low metastatic tumors associate with mutations in GALNTL5, RGPD4, DNAJC5B, DEFA3, IRF2, COL2A1, VEGFC, WFDC11, PROL1, CD80, ADAM18, APOBEC3H and CSHL1; whereas, CYT-high metastatic tumors with PPP6C, DSG1, DEFB112, SPATA16, STARD6, RARRES2 PDE1A, IQCF3, KLK8, CDKN2A, SPANXN5, C16orf90, CD48, TXNDC3 and RAC1.
Although many of these mutations have been previously characterized in skin melanoma (e.g., BRAF, NRAS, TP53, PPP6C, CDKN2A, RAC1, etc.) [10, 80], this is the first report that appreciates them in association with tumors of different cytolytic activity.
We also provide evidence that different genomic structural variations are implicated in the progression of distinct cytolytic subgroups in skin melanoma. We show that CYT-low skin melanomas have a higher number of amplifications and deletions in their genome. These, include mainly recurrent NOTCH2 amplifications, and non-silent mutations and/or deletions in PD-L1, PD-L2, CDKN2A/B, PAX5, ETS1, BCL7, RAD51, JAK2, APAF1, FOXO3, CTNNA1 and IGF1, among others. These genomic structural variations have a profound effect on immune activation in these tumors.
Other distinct chromosomal aberrations included JAK2, KCNH2, MED21 and R3HDM2 amplifications or NRAS, RHOC, CDKN2B, ATM, CHEK1 and ETS1 deletions, and were associated with CYT-high tumors.
Many of them were previously appreciated in the genome of skin melanomas (e.g., MITF, NOTCH2, PD-L1 and JAK2 amplifications, or CDKN2A deletions [10]); however, the CNAs in genomic locations harboring PAX5, ETS1, BCL7, RAD51, APAF1, FOXO3 and CTNNA1 are novel.
Following these observations, we detected marked differences in the chromothriptic pattern between CYT-high and CYT-low tumors. These regarded the prevalence of chromothriptic events, the number of chromothriptic chromosomes per tumor and the cancer genes that they harbor, or the involvement of telomeres and centromeres in the chromothriptic regions. Such differences suggest that distinct mechanisms could give rise to chromothripsis across the two cytolytic subgroups in skin melanoma.
In addition, we highlight several targetable immune-related genes, other than PD-1, PD-L1/2 and CTLA-4, which may coordinately contribute to immunosuppression. IDO1 and LAG3 overexpression, which was evident in CYT-high skin melanomas, provides an alternative immunosuppressive barrier that is needed to hamper the antitumor activity of CTL and NK cells, in these tumors. This is something that was initially proposed in skin melanoma cell lines [110], and later corroborated at the human level, by us and others [14, 36]. We suggest that combined immune inhibition of these markers along with the PD-1/PD-L1 and/or the CTLA-4 axis, could provide a better therapeutic outcome for these patients [19].
Together with the synergistic effect of GZMA and PRF1 in patient survival, we provide evidence that CYT-high skin melanoma patients who received CTLA-4 or PD1/PD-L1/PD-L2 checkpoint blockade therapy (or a combination of both), have a markedly higher immunophenoscore and consequently, an expected clinical benefit, compared to CYT-low patients who did not receive checkpoint inhibition.
In the aggregate, we highlight new links between certain expressional or genomic changes and the activation of an antitumor immune response in skin melanoma. Our findings also corroborate that combination therapies with PD-1, PD-L1/2 and/or CTLA-4 blockade, should be able to overcome resistance and broaden the clinical benefit for these patients.
TCGA, The Cancer Genome Atlas
SKCM, skin melanoma
CYT, immune cytolytic activity
GZMA, granzyme A
PRF1, perforin-1
CTLA-4, cytotoxic T-lymphocyte-associated protein 4
PD-L1, Programmed death-ligand 1
PD-1, Programmed cell death protein 1
NK cells, natural killer cells
CTL, Cytotoxic T cells
SMGs, significantly mutated genes
VAF, variant allele frequency
MATH, mutant-allele tumor heterogeneity
MSigDB, Molecular Signatures Database
CNA, copy number alteration
CDN, classically defined neoepitope
ADN, alternatively defined neoepitope
Ethics approval and consent to participate: All data were retrieved from TCGA’s open access, therefore, no informed consent was necessary.
Consent for publication: Not applicable.
Availability of data and materials: All data generated or analyzed during this study are included in this published article [and its supplementary information files].
Competing interests: The authors declare that they have no competing interests.
Funding: No funding.
Author Contributions: AZ and CR developed methodology. CR, IGS and AZ acquired data, analyzed and interpreted them. AZ designed and supervised the study. AZ also wrote the paper.
Acknowledgements: We are grateful to Mr. Ioannis Apostolides for providing technical assistance with the antigen.garnish R package, as well as the contributors of the TCGA-SKCM dataset and the corresponding skin melanoma patients whose genetic and clinical data were analyzed.
- Bittner M, et al. Molecular classification of cutaneous malignant melanoma by gene expression profiling. Nature. 2000;406(6795):536–40.
- Brash DE, et al. A role for sunlight in skin cancer: UV-induced p53 mutations in squamous cell carcinoma. Proc Natl Acad Sci U S A. 1991;88(22):10124–8.
- Davies H, et al. Mutations of the BRAF gene in human cancer. Nature. 2002;417(6892):949–54.
- Huang FW, et al. Highly recurrent TERT promoter mutations in human melanoma. Science. 2013;339(6122):957–9.
- Leiter U, Garbe C. Epidemiology of melanoma and nonmelanoma skin cancer–the role of sunlight. Adv Exp Med Biol. 2008;624:89–103.
- Lin WM, et al. Modeling genomic diversity and tumor dependency in malignant melanoma. Cancer Res. 2008;68(3):664–73.
- Ruiz R, et al., The RhoJ-BAD signaling network: An Achilles' heel for BRAF mutant melanomas. 2017. 13(7): p. e1006913.
- Winsey SL, et al. A variant within the DNA repair gene XRCC3 is associated with the development of melanoma skin cancer. Cancer Res. 2000;60(20):5612–6.
- León-Letelier RA, Bonifaz LC, Fuentes-Pananá EM. OMIC signatures to understand cancer immunosurveillance and immunoediting: Melanoma and immune cells interplay in immunotherapy. J Leukoc Biol. 2019;105(5):915–33.
- Genomic Classification of Cutaneous Melanoma. Cell, 2015. 161(7): p. 1681–96.
- Danhier F, Feron O, Préat V. To exploit the tumor microenvironment: Passive and active tumor targeting of nanocarriers for anti-cancer drug delivery. J Control Release. 2010;148(2):135–46.
- Fu Q, et al. Prognostic value of tumor-infiltrating lymphocytes in melanoma: a systematic review and meta-analysis. Oncoimmunology. 2019;8(7):1593806.
- Kluger HM, et al. Characterization of PD-L1 Expression and Associated T-cell Infiltrates in Metastatic Melanoma Samples from Variable Anatomic Sites. Clin Cancer Res. 2015;21(13):3052–60.
- Rooney MS, et al. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 2015;160(1–2):48–61.
- Mahoney KM, Freeman GJ, McDermott DF. The Next Immune-Checkpoint Inhibitors: PD-1/PD-L1 Blockade in Melanoma. Clin Ther. 2015;37(4):764–82.
- Hodi FS, et al. Improved survival with ipilimumab in patients with metastatic melanoma. N Engl J Med. 2010;363(8):711–23.
- Luke JJ, et al. Targeted agents or immuno-oncology therapies as first-line therapy for BRAF-mutated metastatic melanoma: a real-world study. Future Oncol. 2019;15(25):2933–42.
- Topalian SL, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med. 2012;366(26):2443–54.
- Christodoulou MI, Zaravinos A. New Clinical Approaches and Emerging Evidence on Immune-Checkpoint Inhibitors as Anti-Cancer Therapeutics: CTLA-4 and PD-1 Pathways and Beyond. Crit Rev Immunol. 2019;39(5):379–408.
- Beatty GL, et al. First-in-Human Phase I Study of the Oral Inhibitor of Indoleamine 2,3-Dioxygenase-1 Epacadostat (INCB024360) in Patients with Advanced Solid Malignancies. Clin Cancer Res. 2017;23(13):3269–76.
- Rasku MA, et al. Transient T cell depletion causes regression of melanoma metastases. J Transl Med. 2008;6:12.
- Rech AJ, Vonderheide RH. Clinical use of anti-CD25 antibody daclizumab to enhance immune responses to tumor antigen vaccination by targeting regulatory T cells. Ann N Y Acad Sci. 2009;1174:99–106.
- Melero I, et al. Agonist antibodies to TNFR molecules that costimulate T and NK cells. Clin Cancer Res. 2013;19(5):1044–53.
- Moran AE, Kovacsovics-Bankowski M, Weinberg AD. The TNFRs OX40, 4-1BB, and CD40 as targets for cancer immunotherapy. Curr Opin Immunol. 2013;25(2):230–7.
- Nirschl CJ, Drake CG. Molecular pathways: coexpression of immune checkpoint molecules: signaling pathways and implications for cancer immunotherapy. Clin Cancer Res. 2013;19(18):4917–24.
- Nicholas C, Lesinski GB. Immunomodulatory cytokines as therapeutic agents for melanoma. Immunotherapy. 2011;3(5):673–90.
- Ngiow SF, et al. Anti-TIM3 antibody promotes T cell IFN-γ-mediated antitumor immunity and suppresses established tumors. Cancer Res. 2011;71(10):3540–51.
- John LB, et al. Anti-PD-1 antibody therapy potently enhances the eradication of established tumors by gene-modified T cells. Clin Cancer Res. 2013;19(20):5636–46.
- John LB, et al. Oncolytic virus and anti-4-1BB combination therapy elicits strong antitumor immunity against established cancer. Cancer Res. 2012;72(7):1651–60.
- Woo SR, et al. Immune inhibitory molecules LAG-3 and PD-1 synergistically regulate T-cell function to promote tumoral immune escape. Cancer Res. 2012;72(4):917–27.
- Takeda K, Okumura K, Smyth MJ. Combination antibody-based cancer immunotherapy. Cancer Sci. 2007;98(9):1297–302.
- Christofi T, et al., Current Perspectives in Cancer Immunotherapy. 2019. 11(10).
- Gide TN, et al. Primary and Acquired Resistance to Immune Checkpoint Inhibitors in Metastatic Melanoma. Clin Cancer Res. 2018;24(6):1260–70.
- Mermel CH, et al. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12(4):R41.
- Zaravinos A, et al. Cytolytic activity correlates with the mutational burden and deregulated expression of immune checkpoints in colorectal cancer. J Exp Clin Cancer Res. 2019;38(1):364.
- Roufas C, et al. The Expression and Prognostic Impact of Immune Cytolytic Activity-Related Markers in Human Malignancies: A Comprehensive Meta-analysis. Front Oncol. 2018;8:27.
- Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
- Law CW, et al. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol. 2014;15(2):R29.
- Law CW, Alhamdoosh M. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. 2016. 5.
- Uhlén M, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.
- Thul PJ, Åkesson L. A subcellular map of the human proteome. 2017. 356(6340).
- Uhlen M, Zhang C. A pathology atlas of the human cancer transcriptome. 2017. 357(6352).
- Finotello F, et al., Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq data. 2019. 11(1): p. 34.
- Newman AM, Liu CL, Green MR. Robust enumeration of cell subsets from tissue expression profiles. 2015;12(5):453–7.
- Carter SL, et al. Absolute quantification of somatic DNA alterations in human cancer. Nat Biotechnol. 2012;30(5):413–21.
- Lawrence MS, et al. Discovery and saturation analysis of cancer genes across 21 tumour types. Nature. 2014;505(7484):495–501.
- Lawrence MS, et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature. 2013;499(7457):214–8.
- Beroukhim R, et al. The landscape of somatic copy-number alteration across human cancers. Nature. 2010;463(7283):899–905.
- Yang J, et al. CTLPScanner: a web server for chromothripsis-like pattern detection. Nucleic Acids Res. 2016;44(W1):W252-8.
- Olshen AB, et al. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics. 2004;5(4):557–72.
- Forbes SA, et al. COSMIC: exploring the world's knowledge of somatic mutations in human cancer. Nucleic Acids Res. 2015;43(Database issue):D805-11.
- Campbell PJ. Telomeres and cancer: from crisis to stability to crisis to stability. Cell. 2012;148(4):633–5.
- Maciejowski J, de Lange T. Telomeres in cancer: tumour suppression and genome instability. Nat Rev Mol Cell Biol. 2017;18(3):175–86.
- Maciejowski J, et al. Chromothripsis and Kataegis Induced by Telomere Crisis. Cell. 2015;163(7):1641–54.
- Carreno BM, et al. Cancer immunotherapy. A dendritic cell vaccine increases the breadth and diversity of melanoma neoantigen-specific T cells. Science. 2015;348(6236):803–8.
- Gros A, et al. Prospective identification of neoantigen-specific lymphocytes in the peripheral blood of melanoma patients. Nat Med. 2016;22(4):433–8.
- Charoentong P, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017;18(1):248–62.
- Alexandrov LB, et al. Signatures of mutational processes in human cancer. Nature. 2013;500(7463):415–21.
- Díaz-Gay M, et al., Mutational Signatures in Cancer (MuSiCa): a web application to implement mutational signatures analysis in cancer samples. 2018. 19(1): p. 224.
- Nik-Zainal S, et al. Mutational processes molding the genomes of 21 breast cancers. Cell. 2012;149(5):979–93.
- Lada AG, et al., AID/APOBEC cytosine deaminase induces genome-wide kataegis. Biol Direct, 2012. 7: p. 47; discussion 47.
- Amelio I, et al., SynTarget: an online tool to test the synergetic effect of genes on survival outcome in cancer. 2016. 23(5): p. 912.
- Bolli N, et al. Heterogeneity of genomic evolution and mutational profiles in multiple myeloma. Nat Commun. 2014;5:2997.
- Chen X, et al. Recurrent somatic structural variations contribute to tumorigenesis in pediatric osteosarcoma. Cell Rep. 2014;7(1):104–12.
- Sakofsky CJ, et al. Break-induced replication is a source of mutation clusters underlying kataegis. Cell Rep. 2014;7(5):1640–8.
- Taylor BJ, et al. DNA deaminases induce break-associated mutation showers with implication of APOBEC3B and 3A in breast cancer kataegis. Elife. 2013;2:e00534.
- Boichard A, Tsigelny IF, Kurzrock R. High expression of PD-1 ligands is associated with kataegis mutational signature and APOBEC3 alterations. Oncoimmunology. 2017;6(3):e1284719.
- Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.
- Ritchie ME, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
- Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.
- Ikeda H, Old LJ, Schreiber RD. The roles of IFN gamma in protection against tumor development and cancer immunoediting. Cytokine Growth Factor Rev. 2002;13(2):95–109.
- Ayers M, et al. IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. J Clin Invest. 2017;127(8):2930–40.
- Han SB, et al. Pectic polysaccharide isolated from Angelica gigas Nakai inhibits melanoma cell metastasis and growth by directly preventing cell adhesion and activating host immune functions. Cancer Lett. 2006;243(2):264–73.
- Kepp O, et al. Extracellular nucleosides nucleotides as immunomodulators. 2017;280(1):83–92.
- Waugh KA, Leach SM, Slansky JE. Tolerance of Tumor-Specific T cells in Melanoma Metastases. J Clin Cell Immunol, 2016. 7(2).
- Bindea G, et al. Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity. 2013;39(4):782–95.
- Girardi M, et al. Regulation of cutaneous malignancy by gammadelta T cells. Science. 2001;294(5542):605–9.
- Ueno H. T follicular helper cells in human autoimmunity. Curr Opin Immunol. 2016;43:24–31.
- Fujimura T, Kambayashi Y, Aiba S. Crosstalk between regulatory T cells (Tregs) and myeloid derived suppressor cells (MDSCs) during melanoma growth. Oncoimmunology. 2012;1(8):1433–4.
- Hodis E, et al. A landscape of driver mutations in melanoma. Cell. 2012;150(2):251–63.
- Dabas N, et al. Diagnostic role of chromosomal instability in melanoma. J Skin Cancer. 2012;2012:914267.
- Hirsch D, et al. Chromothripsis and focal copy number alterations determine poor outcome in malignant melanoma. Cancer Res. 2013;73(5):1454–60.
- Shain AH, et al. The Genetic Evolution of Melanoma from Precursor Lesions. N Engl J Med. 2015;373(20):1926–36.
- Stephens PJ, et al. Massive genomic rearrangement acquired in a single catastrophic event during cancer development. Cell. 2011;144(1):27–40.
- Brás A, Rodrigues AS, Rueff J. Copy number variations and constitutional chromothripsis (Review). Biomed Rep. 2020;13(3):11.
- Maciejowski J, et al., APOBEC3-dependent kataegis and TREX1-driven chromothripsis during telomere crisis. 2020.
- Ernst A, et al. Telomere dysfunction and chromothripsis. Int J Cancer. 2016;138(12):2905–14.
- Zhang CZ, et al. Chromothripsis from DNA damage in micronuclei. Nature. 2015;522(7555):179–84.
- Ly P, et al., Selective Y centromere inactivation triggers chromosome shattering in micronuclei and repair by non-homologous end joining. 2017. 19(1): p. 68–75.
- Schumacher TN, Schreiber RD. Neoantigens in cancer immunotherapy. Science. 2015;348(6230):69–74.
- McGranahan N, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science. 2016;351(6280):1463–9.
- Yarchoan M, et al. Targeting neoantigens to augment antitumour immunity. Nat Rev Cancer. 2017;17(4):209–22.
- Rizvi NA, et al. Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer. Science. 2015;348(6230):124–8.
- Spencer KR, et al. Biomarkers for Immunotherapy: Current Developments and Challenges. Am Soc Clin Oncol Educ Book. 2016;35:e493–503.
- Van Allen EM, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350(6257):207–11.
- Hugo W, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell. 2017;168(3):542.
- Dunn GP, Old LJ, Schreiber RD. The immunobiology of cancer immunosurveillance and immunoediting. Immunity. 2004;21(2):137–48.
- Gajewski TF, Louahed J, Brichard VG. Gene signature in melanoma associated with clinical activity: a potential clue to unlock cancer immunotherapy. Cancer J. 2010;16(4):399–403.
- Harlin H, et al. Chemokine expression in melanoma metastases associated with CD8 + T-cell recruitment. Cancer Res. 2009;69(7):3077–85.
- Spranger S, et al. Up-regulation of PD-L1, IDO, and T(regs) in the melanoma tumor microenvironment is driven by CD8(+) T cells. Sci Transl Med. 2013;5(200):200ra116.
- Fridman WH, et al. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12(4):298–306.
- Denkert C, et al. Tumor-associated lymphocytes as an independent predictor of response to neoadjuvant chemotherapy in breast cancer. J Clin Oncol. 2010;28(1):105–13.
- Tumeh PC, et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature. 2014;515(7528):568–71.
- Selitsky SR, et al. Prognostic value of B cells in cutaneous melanoma. Genome Med. 2019;11(1):36.
- Falleni M, et al. M1 and M2 macrophages' clinicopathological significance in cutaneous melanoma. Melanoma Res. 2017;27(3):200–10.
- Ladányi A, et al. Density of DC-LAMP(+) mature dendritic cells in combination with activated T lymphocytes infiltrating primary cutaneous melanoma is a strong independent prognostic factor. Cancer Immunol Immunother. 2007;56(9):1459–69.
- Pietra G, et al. Natural killer cells kill human melanoma cells with characteristics of cancer stem cells. Int Immunol. 2009;21(7):793–801.
- Balli D, et al. Immune Cytolytic Activity Stratifies Molecular Subsets of Human Pancreatic Cancer. Clin Cancer Res. 2017;23(12):3129–38.
- Morrison C, et al. Predicting response to checkpoint inhibitors in melanoma beyond PD-L1 and mutational burden. J Immunother Cancer. 2018;6(1):32.
- Pietra G, et al. Melanoma cells inhibit natural killer cell function by modulating the expression of activating receptors and cytolytic activity. Cancer Res. 2012;72(6):1407–15.