Analysis of DNA methylation patterns in the tumor immune microenvironment of metastatic melanoma.

The presence of immune cells in the tumor microenvironment has been associated with response to immunotherapies across several cancer types, including melanoma. Despite its therapeutic relevance, characterization of the melanoma immune microenvironments remains insufficiently explored. To distinguish the immune microenvironment in a cohort of 180 metastatic melanoma clinical specimens, we developed a method using promoter CpG methylation of immune cell type‐specific genes extracted from genome‐wide methylation arrays. Unsupervised clustering identified three immune methylation clusters with varying levels of immune CpG methylation that are related to patient survival. Matching protein and gene expression data further corroborated the identified epigenetic characterization. Exploration of the possible immune exclusion mechanisms at play revealed likely dependency on MITF protein level and PTEN loss‐of‐function events for melanomas unresponsive to immunotherapies (immune‐low). To understand whether melanoma tumors resemble other solid tumors in terms of immune methylation characteristics, we explored 15 different solid tumor cohorts from TCGA. Low‐dimensional projection based on immune cell type‐specific methylation revealed grouping of the solid tumors in line with melanoma immune methylation clusters rather than tumor types. Association of survival outcome with immune cell type‐specific methylation differed across tumor and cell types. However, in melanomas immune cell type‐specific methylation was associated with inferior patient survival. Exploration of the immune methylation patterns in a pan‐cancer context suggested that specific immune microenvironments might occur across the cancer spectrum. Together, our findings underscore the existence of diverse immune microenvironments, which may be informative for future immunotherapeutic applications.

The presence of immune cells in the tumor microenvironment has been associated with response to immunotherapies across several cancer types, including melanoma. Despite its therapeutic relevance, characterization of the melanoma immune microenvironments remains insufficiently explored. To distinguish the immune microenvironment in a cohort of 180 metastatic melanoma clinical specimens, we developed a method using promoter CpG methylation of immune cell type-specific genes extracted from genome-wide methylation arrays. Unsupervised clustering identified three immune methylation clusters with varying levels of immune CpG methylation that are related to patient survival. Matching protein and gene expression data further corroborated the identified epigenetic characterization. Exploration of the possible immune exclusion mechanisms at play revealed likely dependency on MITF protein level and PTEN loss-of-function events for melanomas unresponsive to immunotherapies (immune-low). To understand whether melanoma tumors resemble other solid tumors in terms of immune methylation characteristics, we explored 15 different solid tumor cohorts from TCGA. Low-dimensional projection based on immune cell type-specific methylation revealed grouping of the solid tumors in line with melanoma immune methylation clusters rather than tumor types. Association of survival outcome with immune cell typespecific methylation differed across tumor and cell types. However, in melanomas immune cell type-specific methylation was associated with inferior patient survival. Exploration of the immune methylation patterns in a pan-cancer context suggested that specific immune microenvironments might occur across the cancer spectrum. Together, our findings underscore the existence of diverse immune microenvironments, which may be informative for future immunotherapeutic applications.
The tumor microenvironment is a complex mixture of malignant and nonmalignant cells that includes immune cells. Clearly, the microenvironment of melanoma tumors, especially immune cells, plays an important role in determining progression (Passarelli et al., 2017;Tucci et al., 2018;Villanueva and Herlyn, 2008). Single-cell RNA sequencing of melanoma tumors has confirmed that a complex mixture of immune T-cell subtypes exists within a tumor (Jerby-Arnon et al., 2018;Sade-Feldman et al., 2018). In addition, bulk transcriptomic analyses have been utilized to explore the immune microenvironment from melanoma tumors, demonstrating that tumors with increased expression of genes involved in a diverse range of immune systems have increased survival (Angelova et al., 2015;Bindea et al., 2013). Current methods to decipher the immune landscape in tumors are based on the transcriptome (Angelova et al., 2015;Bindea et al., 2013). However, global DNA methylation patterns are largely associated with cellular lineage and show higher distinction between cellular lineages than mRNA expression in blood and skin lineages (Bock et al., 2012). Moreover, DNA methylation has been demonstrated to resolve cell of origin of peripheral blood cells (Houseman et al., 2012) and cell-free DNA (Moss et al., 2018), and was introduced as a complementary approach to classify central nervous system (CNS) tumors (Capper et al., 2018). Nevertheless, DNA methylation has not been extensively exploited to determine the immune lineages occurring in the microenvironment of melanoma. Therefore, we created DNA methylation signatures that reflect a wide range of tumor-associated immune cell subsets. These signatures grouped melanoma tumors into three clusters with distinct clinical and molecular properties that were re-identified in independent data. Further molecular exploration revealed several immune exclusion mechanisms activated in tumors. Moreover, an immune-rich microenvironment consisting of a broad range of immune cells was observed across several cancer types. In conclusion, DNA immune methylation signatures harbor significant clinical and biological information.

Patient sample information
The study was approved by the Regional Ethics Committee at Lund University (Dnr. 191/2007 and101/ 2013). The experiments were undertaken with the understanding and written consent of each subject. The study methodologies conformed to the standards set by the Declaration of Helsinki. The sample cohort consisted of 214 melanoma tumors and was obtained at the Department of Surgery at Sk ane University Hospital. This is a historic cohort collected before the era of ICIs and targeted therapies.

Sample preparation and bisulfite conversion
DNA extraction from patient samples was done as described previously (Harbst et al., 2014) and bisulfiteconverted using EZ DNA Methylation Kit (Zymo Research Europe, Freiburg, Germany) according to the manufacturer's instructions and analyzed using Illumina Infinium MethylationEPIC BeadChip array. Methylation data are available in Gene Expression Omnibus with accession number GSE144487.

Immune cell type-specific CpG set for metastatic melanoma
We selected two immune cell type-associated gene signatures (Angelova et al., 2015;Tirosh et al., 2016)  identify immune cell types in melanoma tumors. Next, we combined both gene signatures, and after removing conflicting genes (genes belonging to different cell types in each signature), we obtained a total of 920 genes associated with 30 immune cell subsets. Next, we selected corresponding promoter CpGs (< 1500 bp from TSS) of these genes, since their association with the transcriptional regulation of a gene is well known and thus more likely to be lineage-specific. After mapping the gene signatures with the list of CpG probes shared between 450K and EPIC platforms, we found 4228 promoter CpGs belonging to 744 genes and associated with 30 immune cell subsets.
The CpGs were then further filtered using the following criteria: 1 Selecting the most differentially methylated CpGs across reference immune cells: We dichotomized the b-values into robust methylation bins, as unmethylated (b < 0.3) and methylated (b ≥ 0.3). We then selected CpGs that have significantly different proportions of methylated and unmethylated signals among the reference immune cells using Fisher's exact test and an FDR < 0.01. 2 High methylation in nonimmune cells: To ensure that any methylation difference we observe is likely coming from immune cells and not from other cells present in the microenvironment, we further shortlisted CpGs with a high percentage and level of methylation (b > 0.7 in > 98% of samples) among nonimmune normal cells and melanoma cell lines. 3 Forming gene-CpG pairs: We wanted to ensure that any single gene is not over-represented through the presence of multiple CpGs. We therefore selected the most significant CpG for each gene from Fisher's exact test in step 1.
The selection processes resulted in 67 gene-CpG pairs belonging to 21 immune cell populations.

Immune cell type-specific CpG set for nonmelanoma TCGA pan-cancer cohorts
The CpG selection procedure was identical to the process we followed for metastatic melanoma (MM) tumor cohorts except step 2. At step 2, we filtered the immune CpGs against methylation profiles of matching tumor cell lines from Genomics of Drug Sensitivity (GDSC) database and selected CpGs for further analyses if they had shown high level and percentage of methylation in the tumor cell lines (b > 0.7 in > 90% of samples). Here, we had to relax the sample selection criteria since for a number of GDSC tumor types a minority of cell lines were displaying a methylation pattern that deviated from the majority of cell lines of the tumor type.

Immune methylation centroid-based classification
Immune methylation centroid-based classification was performed by correlating sample methylation profiles across centroid CpGs to each cluster centroid (Table S3) and then selecting the cluster that reported highest correlation (s Kendall ). If no cluster displayed correlation ≥ 0.3, then the sample was annotated unclassified.
2.6. Immune cell type-specific methylation score calculation Immune cell type-specific methylation scores were calculated using the matching CpGs from the 67 CpG set for MM cohorts and by taking median methylation value of all CpGs belonging to a specific immune cell type. For non-MM solid tumor cohorts, the process is almost identical to MM but here cohort-specific immune CpG sets were used for score calculation.

PTEN promoter hypermethylation calculation
We selected promoter CpGs for the PTEN gene that are located at the DNase hypersensitivity sites (DHS) of the promoter, as the PTEN promoter contained a complex set of CpGs on the Illumina 450K array. Next, we called PTEN hypermethylation in tumors if more than 10% of the DHS promoter CpGs are hypermethylated (b > 0.7) and median b-value for all PTEN promoter CpGs for the corresponding tumor is above 0.5.

Statistical analyses and calculations for immune cell type methylation and gene expression scores
All statistical and bioinformatics analyses were performed in R. For comparing numerical values, we used Spearman and Kendall correlation. Comparisons between two groups were performed using Mann-Whitney U-test/Wilcoxon rank-sum test, and for more than two groups, we used Kruskal-Wallis test. For survival analyses, we used univariate and multivariate Cox regression.
Both methylation and gene expression scores were calculated by median methylation/expression of CpGs/genes belonging to each immune cell type. We used median-centered gene expression for calculating expression scores in Lund and TCGA cohort analyses; noncentered gene expression scores were used in pancancer analysis of correlation between methylation and gene expression.

DNA methylation-based immune profiling of metastatic melanoma tumors
To explore the DNA methylation immune landscape in MM tumors, we first curated two previously described immune gene expression signatures (Angelova et al., 2015;Tirosh et al., 2016) (Fig. 1A) to identify corresponding gene promoter CpGs. Next, we identified those CpGs that were most differentially methylated across seven normal peripheral bloodderived immune cell subsets (Reinius et al., 2012), while in parallel showing high methylation in nonimmune cell types in the tumor microenvironment (cancer cell lines and normal stromal cells-dermal fibroblasts, epidermal keratinocytes, and melanocytes) . After subsequent selection of the best CpG (CpG showing lowest FDR in Fisher's exact test for differential methylation) for individual genes, we arrived at a final list of 67 CpGs (Table S1). These, we hypothesize, are DNA methylation signatures representative of different immunological cell types (Fig. 1B). Next, we used this CpG set to analyze Illumina EPIC methylation array-based profiles from 180 MM tumors (the Lund cohort, Table S2).

Immune cell type-associated CpGs reveal melanoma subgroups
Consensus clustering of the Lund cohort using the selected immune cell type-specific CpGs identified three immune methylation clusters (Fig. S1A). Cluster 1 showed overall lower methylation in the majority of the immune CpGs, whereas Cluster 3 displayed the opposite, with Cluster 2 being intermediate in appearance ( Fig. 2A). To explore the observed methylation differences further, we constructed a median methylation score for all CpGs associated with the described immune cell types for each tumor. We found that differences in methylation scores between Clusters 1 and 3 were most pronounced, with 19 out of 21 immune cell types being significant (Fig. 2B, Kruskal-Wallis and post hoc Dunn test; multiple testing corrections across cell types for each comparison individually, FDR < 0.001). The biggest difference in methylation values was observed for central memory CD4 + T cells (absolute Db: 1 vs 2: 0.25; 2 vs 3: 0.14; 1 vs 3: 0.39; FDR < 0.001) and effector memory CD8 + T cells (absolute Db: 1 vs 2: 0.3; 2 vs 3: 0.18; 1 vs 3: 0.48; FDR < 0.001). Significant differences between Clusters 1 and 2 were observed mostly in the lymphoid lineage and differences between Clusters 2 and 3 followed a similar pattern of significance to that of Clusters 1 and 3, albeit with lower level of methylation difference (absolute Db-value).

Immune methylation clusters are associated with patient survival
To determine whether the identified methylation patterns have prognostic value, we analyzed for differences in the survival outcome of the melanoma patients in Lund cohort. Immune methylation clusters showed significant difference in survival using both distant metastasis-free survival (DMFS) for regional metastatic cases and melanoma-specific survival for the entire metastatic cohort (DSS) using univariate Cox regression models (Fig. 2C,D, respectively). Cluster 2 showed significantly decreased DMFS [hazard ratio (HR): 2.45; 95% confidence interval (CI): 1.30-4.61; P = 0.005] compared to Cluster 1. Regarding DSS, we found both Cluster 2 (HR: 2.10; 95% CI: 1.22-3.64; P = 0.008) and Cluster 3 (HR: 2.18; 95% CI: 1.33-3.56; P = 0.002) to be of significantly higher risk when compared to Cluster 1. Furthermore, multivariable Cox regression models with adjustments for age and gender (DMFS) and metastasis type (DSS) reported analogous findings for both DMFS (Fig. S1B) and DSS (Fig. S1C). Overall, using DNA methylationbased immune profiling we identified three melanoma subgroups with different patient survival.

Molecular characterization of immune methylation clusters
Next, we set out to explore whether the identified immune methylation clusters were associated with transcriptional differences of the immune cell-specific genes in matched tumors (n = 179). Overall, transcriptional differences followed the observations on DNA methylation levels (Fig. 3A). Indeed, for most selected immune cell CpGs, moderate-to-high negative correlations (À0.2 to À0.5 and below, s Kendall ) were found between methylation and expression of the corresponding genes, as expected for the promoter CpGs involved in the regulation of gene expression (Fig. S1D). Next, we investigated whether the identified methylation clusters were accompanied by differences in immune cell infiltration using hematoxylin and eosin staining for tumor-infiltrating lymphocytes (TILs) along with immunohistochemical staining (IHC) for CD3 and CD8 for T cells, and CD68 and CD163 for myeloid cells. Using matched staining data from 127 of 180 MM tumors, we analyzed association between clusters and staining categories (strong infiltration, localized infiltration, and absent for TILs and T cells; absent/ low infiltration, nontumor infiltration, and tumor infiltration for myeloid cells). In line with the DNA methylation data, Cluster 1 had the highest percentage of samples with strong and localized infiltration of TILs, CD3 + , and CD8 + T cells (P < 0.001; FDR < 0.001) (Fig. 3B). The presence of TILs showed significant association with DSS in the Lund cohort ( Fig. S1E). Regarding CD68 + and CD163 + myeloid cells, Cluster 1 had fewer absent/low-infiltrating samples compared to Cluster 3. Together, these results support that tumor infiltration of lymphocytes and myeloid cells is accurately reflected by DNA methylation.
Recently, the existence and importance of tertiary lymphoid structures (TLS) in the melanoma tumors have been reported along with their close relationship with prognosis and response to immunotherapeutic treatments (Cabrita et al., 2020;Helmink et al., 2020). Here, using the matched TLS information for the tumors (Cabrita et al., 2020), we observed significant enrichment of TLSs in Cluster 1 tumors (P < 0.001; FDR < 0.001) (Fig. 3B), indicating the possibility of a successful antitumor immunity in these tumors.

Immune cell type-associated gene sets from Angelova et al. and Tirosh et al.
Combining

DSS -All metastases
Log-rank P = 0.003 Kaplan-Meier plots of the immune methylation clusters for DMFS (regional lymph node metastases only) and DSS (all metastases except local recurrences), respectively. Numbers at risk and percentages at risk (within brackets) are indicated below each plot. P-values calculated using log-rank test. Maximum time to display on x-axis has been clipped to 10 years. Tumor mutational burden has emerged as an important biomarker for predicting response to ICI treatments (reviewed in Chan et al., 2019). Using mutation data from 1368 cancer genes in 129 of 180 MM tumors from the Lund cohort, we did not observe a relationship between immune methylation clusters and tumor mutational burden or specific melanoma hotspot mutations (Fig. 3C,D, respectively). Additionally, we analyzed association between immune methylation clusters and our previously reported gene expression phenotypes  and observed strong association between high-immune class and low methylation Cluster 1 (Fisher's exact P < 0.001, Table S3).

Survival probability Survival probability
Overall, the immune methylation clusters are independent of mutational patterns and tumor mutational burden but are confirmed by the presence or absence of immune cells on the mRNA and protein level.

Analysis of immune exclusion mechanism associated with immune methylation clusters
Several molecular mechanisms have been proposed to explain why certain tumors are able to escape the immune system (Casey et al., 2016;Gettinger et al., 2017;Le et al., 2017a;Peng et al., 2016;Spranger et al., 2015;further reviewed in Spranger and Gajewski, 2018;Trujillo et al. 2018). Clearly, immune methylation Cluster 3 was characterized by an immune-poor DNA methylation pattern reflected by increased methylation levels of immune cell-specific CpGs ( Fig. 2A). Further, we hypothesized that the differentiation state of the melanoma tumor may determine the ability of Cluster 3 tumors to evade the immune system. Thus, we used matched staining information for 127 of 180 MM tumors for the MITF protein, which is a melanocytespecific transcription factor important in melanocyte differentiation and melanoma cell survival (Carreira et al., 2005;Garraway et al., 2005;McGill et al., 2002;reviewed in Levy et al., 2006;Steingr ımsson et al., 2004). Overall, 25% of the matched stained tumors were negative for MITF protein. Intriguingly, we found fewer MITF-negative tumors in the immune-rich methylation Cluster 1 ( Fig. 4A; chi-square P = 0.06).
As MITF-positive and MITF-negative tumors are known to harbor different transcriptional programs and presumably represent two distinct melanoma states (Hoek et al., 2006), we thus performed subsequent analyses in MITF-positive and MITF-negative tumors separately. In pathway analysis of gene expression data, we found that MITF-positive tumors showed enrichment for c-MYC (MYC) target genes and DNA repair-associated genes in Cluster 3 (FDR < 0.05); however, such observations were not found for the MITF-negative group. Upon stratification by MITF status, we observed a significant difference in MYC expression across clusters for MITF-positive tumors (Fig. 4B, P = 0.001), which was not observed for the combined data (Fig. S2A). Furthermore, we explored the molecular profiles of known immune exclusion-associated genes (beta-catenin pathway, PI3K pathway, TP53, reviewed in Spranger and Gajewski (2018)) across clusters. CTNNB1 showed a significant mRNA expression difference across clusters (Fig. 4C, P < 0.001, FDR < 0.05); however, when stratified by MITF status, characteristics similar to MYC were revealed ( Fig. S2B, P < 0.001). Together, these observations delineate that MYC and CTNNB1 may act as immune exclusion molecules in MITF-positive melanoma tumors. Among PI3K-Akt pathway genes, we observed significant gene expression differences across clusters for PTEN and PIK3R1 (Fig. 4D and Fig. S2C, P = 0.005; FDR < 0.05 and P = 0.03; FDR < 0.05, respectively) irrespective of MITF status. Upon further investigation, we found that PTEN somatic copy number alterations (SCNAs) showed good overlap with PTEN promoter hypermethylation and nonsynonymous mutations for the same tumors. However, PTEN hypermethylation and mutations were found to be mostly exclusive (Fig. 4E, Fig. S2D). Next, we summed up all molecular alteration events of PTEN for each tumor and found significant concordance between PTEN alteration events and corresponding mRNA expression (Fig. S2E). Also, the distribution of tumors harboring any PTEN alteration event varied significantly across clusters (chi-square P < 0.001) with Cluster 3 showing highest enrichment (Fig. 4F). We further explored the association between PTEN mRNA expression and TILs and observed a trend in difference in terms of PTEN expression between tumors with absent and brisk TILs (Fig. S2F). The association of TILs with molecular events in PTEN was more prominent (Fig. 4E), and a significant difference in terms of TILs was observed between tumors with and without PTEN molecular events (chi-square P = 0.001). Association between immune methylation clusters and molecular event status in PTEN was also validated in the TCGA MM cohort (Fig. S2G). Together, we observed a known immune exclusion mechanism in the immune-poor MM tumors displaying a complex pattern with upregulation of PI3K-Akt pathway by blocking PTEN functioning through molecular alterations and upregulation of MYC and CTNNB1 in MITF-positive tumors only.

Validation of immune methylation clusters in external melanoma cohorts
To confirm the immune methylation clusters, we built methylation-based centroids to classify samples belonging to the TCGA MM cohort. For this, we selected CpGs from the immune CpG set (Table S1) which showed a significant methylation difference across Lund clusters (Kruskal-Wallis test, FDR < 0.0001). Median methylation values for these significant CpGs (n = 51) were computed for each cluster to form methylation centroids (Table S4). A minor fraction of samples (2.37%) could not be classified (s Kendall < 0.3). Immune methylation clusters in TCGA cohort displayed methylation profiles similar to corresponding clusters in the Lund cohort (Fig. S3A). When we grouped both the Lund and TCGA cohorts together using common centroid CpGs and tSNE (van der Maaten and Hinton, 2008) (perplexity = 30; h = 0.1), we noticed clear segregation between the immune methylation clusters in the first two dimensions, as expected (Fig. 5A, left). Such segregation was not observed when samples were annotated on the basis of the site of metastasis (metastasis type) or cohort (Fig. 5A middle and right, respectively). Additionally, gene expression scores, tumor mutational burden, and mutational patterns displayed identical characteristics for clusters as observed in the Lund cohort (Fig. 5B,C,  Fig. S3B, respectively). Importantly, Cluster 3 displayed a significant inferior overall survival (OS, HR: 2.29; 95% CI: 1.54-3.40; P < 0.001) compared to Cluster 1 in TCGA data (Fig. 5D) and remained significant in multivariate analyses with adjustment for metastasis type (P < 0.001, Fig. S3C). Additionally,  we analyzed the association between immune methylation clusters in TCGA cohorts and previously reported gene expression classes for the same cohort (Akbani et al., 2015) and observed strong association between Immune class and low methylation Cluster 1 (Fisher's exact P < 0.001, Table S5).
To further explore immune gene methylation in distant metastases, we classified brain, lymph node, and other metastases from distant sites from TCGA skin cutaneous melanoma (SKCM) cohort (n = 68) and brain metastases (n = 30) from another cohort (GSE108576) (Orozco et al., 2018) into the immune methylation clusters using methylation-based centroids. We found that samples in both cohorts were enriched for Cluster 3 irrespective of the metastasis location (Fig. 5E). In TCGA cohort, we did not find any significant association between metastasis location and cluster membership. When both cohorts were analyzed together across the common centroid CpGs and tSNE (perplexity 30; h = 0.1), observations similar to earlier Lund and TCGA joint analysis were made (Fig. 5F). Together, the immune cell methylation characteristics of the Lund cohort were re-identified in external data, with similar prognostic implications, and the immune methylation profiles were largely independent of the site of metastasis.

Implications of the epigenetic immune clusters in a pan-cancer context
Next, we wanted to investigate whether immune methylation signatures have the same prognostic implications in other cancer types as in melanoma, and to which extent immune environments display shared characteristics among solid tumors. At first, we wanted to dissect the role of specific immune cell lineages across cancer types. To this purpose, we created immune cell methylation scores for MM tumors (both Lund and TCGA cohorts) using previously derived melanoma-specific immune CpG set. Similar scores for non-MM TCGA solid tumor cohorts were constructed by identifying immune cell type-specific CpGs individually for each cohort in a process analogous to MM. Here, the immune CpGs were filtered against cancer type-specific CpGs using cell lines from the Genomics of Drug Sensitivity (GDSC) database. The resulting methylation scores were highly anticorrelated to the corresponding gene expression scores (computed in the similar manner as methylation scores, using genes instead of the corresponding CpGs), for most immune cell types and tumor cohorts (Fig. 6A). Next, to explore the relations among the immune microenvironments across tumor types, we constructed a tSNE plot using aforementioned methylation scores for 17 immune cell types. When tumor types were annotated on the basis of their tissue of origin, we observed some tumor cohorts were more localized than others (Fig. 6B, Left). Known immunedeprived tumor types like CNS tumors (lower grade glioma, LGG; and glioblastoma multiforme, GBM) colocalized with a group of similarly immune-poor tumor types from the developmental gastrointestinal tumors (pancreatic adenocarcinoma, PAAD; and liver hepatocellular carcinoma, LIHC), thus indicating that immune microenvironments of these tumor types bear resemblance despite their diverse tissue of origin. Moreover, the immune methylation Clusters 1 and 3 from MM separated well in this immune microenvironment plot (Fig. 6B, Right). Cluster 1 tumors were found to be colocalizing with a group of diverse tumor types, indicating existence of similar type of immunologically active tumors across the cancer spectrum. Conversely, we found Cluster 3 tumors colocalizing with most CNS and endocrine tumors, indicating that the immune microenvironment of these immune-poor melanomas likely resembles other immune-poor tumor types. Next, we analyzed the prognostic implications of immune cell type methylation scores across the cohorts using overall patient survival (OS). For this, we dichotomized the methylation scores into two categories, b > 0.7 as hypermethylated and b ≤ 0.7 as not hypermethylated. We found methylation scores for CD56 + dim NK cells, and activated and effector memory CD8 + T cells, to be significantly prognostic across multiple cancer cohorts including melanoma, with low methylation levels inferring good prognosis in most cancer types (Fig. 6C). Overall, we found low methylation scores of immune cell types to be associated with good outcome for the majority of cell types, particularly for immune cells from the lymphoid lineage, which might reflect the role of the adaptive immune system for an effective response to tumor neoantigens. However, the majority of significant associations with good outcome (adjusted P < 0.1) occurred in only two tumor types, that is, MM and head and neck squamous cell carcinoma (HNSC). Collectively, methylation-based immune cell signatures inform on prognosis in several solid tumor cohorts, particularly in melanoma.

Discussion
In this study, we report that immune methylation-specific signatures harbor important prognostic information and simultaneously inform us about the different tumor immune microenvironments across cancer landscape with a particular focus on MM. DNA methylation has established its role as a major epigenetic driver in cancer   progression and development. However, its contribution to define the characteristics of tumor microenvironment remains poorly understood. It has been pointed out recently that DNA hypomethylation promotes immune evasion in the corresponding tumors (Jung et al., 2019). Also, DNA methylation pattern predicting response to ICI treatment for non-small-cell lung cancer has been unveiled (Duruisseaux et al., 2018). Thus, further emphasis on DNA methylation for the distinction of different tumor microenvironments is warranted. Our analyses with an immune cell type-specific CpG set have unraveled diversity in the immune landscape of MM tumors in terms of DNA methylation. We identified three immune methylation clusters that were significantly associated with patient survival and further supported by transcriptomic and immunostaining data. Together, these analyses have hinted at the inverse relationship between immune cell type-specific methylation and enrichment of the corresponding immune cell types in the tumor microenvironment.
Tumor mutational burden has often been associated with positive clinical response to ICI treatments for several tumor types including melanoma (Goodman et al., 2017;Van Allen et al., 2015, reviewed in Chan et al., 2019. However, we did not observe a clear association of immune methylation clusters with the tumor mutational burden, in line with a previous report (Spranger et al., 2016). This might be due to a universally high immunogenicity of melanoma cells, where immune evasion, as observed in Cluster 3, that is, immune methylation high cluster, can rather be achieved by transcriptional changes (Jerby-Arnon et al., 2018). Therefore, in this study we additionally explored the possible immune exclusion mechanisms (reviewed in Spranger and Gajewski, 2018) at play for the Cluster 3 tumors in our cohort. Upregulation of PI3K-Akt pathway and subsequent downregulation of PTEN are a major immune evasion mechanism that has been observed across cancers including melanoma (reviewed in Dituri et al., 2011;Spranger and Gajewski, 2018). However, somatic mutations in PTEN gene have found to be infrequent in melanomas (reviewed in Aguissa-Tour e and Li, 2012); hence, alternative molecular mechanisms such as SCNAs and promoter hypermethylation have been proposed as an alternative mechanism of blocking PTEN functions (Roh et al., 2016;Stahl et al., 2003). Our observations for Cluster 3 tumors further supported the hypothesis that immune-poor melanomas likely achieve immune evasion through blocking of PTEN functioning through promoter hypermethylation and SCNAs and often in conjunction with one another. The role of MITF as a major transcription factor regulating melanoma progression and development is well established (Carreira et al., 2005;Garraway et al., 2005;McGill et al., 2002;reviewed in Levy et al., 2006;Steingr ımsson et al., 2004). MITF-low melanomas have been shown to be resistant to multiple targeted treatments (M€ uller et al., 2014). Studies have also shown that decreased expression of melanocyte differentiation antigens through downregulation of MITF could possibly trigger immune evasion in melanomas (Kono et al., 2006;Woods et al., 2014). Nonetheless, exact role of MITF in the immune exclusion mechanism is still unclear. Our analyses with immune-poor melanomas stratified by MITF immunostaining suggested transcriptional upregulation of b-catenin-dependent canonical Wnt signaling pathway along with the upregulation of c-Myc in the MITF-positive group. Immune evasion through upregulation of canonical Wnt signaling pathway is common in melanomas (Spranger et al., 2015), and the role of c-Myc in the immune exclusion process has previously been described (reviewed in . Since b-catenin-induced melanomas LGG, lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma of lung; PAAD, pancreatic adenocarcinoma; STAD, stomach adenocarcinoma; THCA, thyroid carcinoma). MM cohorts are highlighted in red. Methylation scores for MM cohorts were calculated from the melanoma-specific 67 immune CpG set, whereas for other tumor types, they were calculated using individual tissue type-specific immune CpG sets in the similar fashion to MM. (B) tSNE plot of pan-cancer cohorts based on the previously mentioned methylation scores of common 17 immune cell types; plot is annotated by tissue of origin (left) and immune methylation cluster membership for MM (right). (C) Balloon plot of hazard ratios (HR) for dichotomized immune methylation scores (b > 0.7: hypermethylated, b ≤ 0.7: not hypermethylated) from univariate Cox regression models for each immune cell type and cancer cohorts. MM cohorts are highlighted in red. Balloon colors are indicative of direction of hazard ratios, and sizes indicate value category of HR (large: HR ≥ 1.5 or HR ≤ 0.66; small: HR < 1.5 or HR > 0.66). Additional black circles indicate significance after correction for multiple testing (Benjamini-Hochberg corrected P-value < 0.1). Gray circles indicate situations where HR calculations were not possible, and missing circles indicate cell types where no CpG survived filtering.  Lund SKCM  TCGA SKCM  TCGA THCA  TCGA CRC  TCGA LUAD  TCGA LUSC  TCGA ESCA  TCGA HNSC  TCGA GBM  TCGA KIRC  TCGA LGG  TCGA LIHC  TCGA BLCA  TCGA STAD  TCGA BRCA  TCGA MESO  TCGA   immunotherapeutic treatments (Cristescu et al., 2018;Van Allen et al., 2015), hence these results suggest further studies to clarify the role of PTEN and MITF levels in immune exclusion. Notably, our analyses on distant MMs did not reveal a major difference in the immune microenvironment characterization among tumors from different locations, most notably brain. This is further supported by clinical trials showing efficacy of ICI for intracranial metastases concordant with the extracranial ones and systemic response (Kluger et al., 2019;Margolin et al., 2012;Tawbi et al., 2018). A recent study using mouse models suggested such efficacy is likely due to infiltration of T cells especially CD44 + CD62L À effector memory cells from the extracranial sites (Taggart et al., 2018).

T C G A -S K C M T C G A -B L C A T C G A -B R C A T C G A -C R C T C G A -E S C A T C G A -G B M T C G A -H N S C T C G A -K I R C T C G A -L G G T C G A -L I H C T C G A -L U A D T C G A -L U S C T C G A -M E S O T C G A -P A A D T C G A -S T A D T C G A -T H C A
Our analyses on a pan-cancer cohort revealed similarities in the immune microenvironments of tumors from diverse tumor types in terms of immune methylation. Generally, our findings are compatible with the concept of inflammatory and noninflammatory tumor environments (Chakravarthy et al., 2018); however, our pan-cancer analyses suggest that immune environments carry additional complexity beyond this concept, as has been demonstrated before (Thorsson et al., 2018). Using cohort-specific immune methylation scores, we observed grouping of MM tumors belonging to the low immune methylation cluster along with tumors from lung and gastrointestinal (GI) tract. Such observation indicates that immune-rich microenvironments from melanoma have similarities to microenvironment of lung and GI tumors. In this regard, it is interesting that ICI treatments have garnered significant attention for non-small-cell lung cancer (NSCLC), which is among the cancer types where such an immune-rich microenvironment likely occurs (reviewed in Moon et al., 2017). Therefore, in addition to a high mutational burden, such as induced by microsatellite instability (MSI) (Le et al., 2015), an immune-rich environment may confer improved outcome to immunotherapy agents. Among the different immune cell types, the prognostic value of effector memory T cells could also be established in several tumor types of the pan-cancer cohort. Memory T cells have received attention recently, as their presence in the tumor microenvironment was associated with clinical response to ICI treatment in melanoma (Gide et al., 2019;Sade-Feldman et al., 2018).

Conclusion
In summary, we explored the immune microenvironment in MM tumors from a DNA methylation perspective. DNA methylation offers several advantages over transcriptomic characterization of mixed cellular environments, primarily due to its higher frequency in terms of differences between cellular lineages compared to gene expression (Bock et al., 2012). Additionally, the biologically limited range of methylation b-values, being restricted from 0 (unmethylated) to 1 (methylated), facilitates the comparison across experiments and laboratories, whereas transcriptomic data naturally have a higher dynamic range. Also, greater stability of DNA compared to RNA in formalin-fixed, paraffin-embedded (FFPE) samples (Okello et al., 2010) makes it easier to work with DNA methylation in archival, historical cohorts.
Our use of a largely untreated historic cohort ensured that the immune cell type methylation patterns are not biased by treatment. Also, the extension of our study to other solid tumors revealed an immune-rich environment in several other cancer types. The prognostic implications of some innate and adaptive immune cell types were re-identified in other cancer types. These results point toward a more tailored approach to immunotherapy based on tumor immune microenvironment. fan Paulsson's Foundation, and the governmental funding for healthcare research (ALF). This study was also supported by the Intramural Research Program (IRP) of the Division of Cancer Epidemiology and Genetics, National Cancer Institute, US National Institutes of Health (NIH); the contents of the publication do not necessarily reflect the views of policies of the US Department of Health and Human Services, nor does mention of trade names, commercial products, or organizations imply endorsement by the US government. This work utilized the computational resources of the NIH high-performance computational capabilities Biowulf cluster (http://hpc.nih.gov). We acknowledge contributions to measurement of methylation from The National Cancer Institute Cancer Genomics Research Laboratory (CGR). This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 641458. MR is financially supported by the Knut and Alice Wallenberg Foundation as part of the National Bioinformatics Infrastructure Sweden at Sci-LifeLab. We also wish to thank the members of the Swedish Bioinformatics Advisory Program and MEL-GEN consortium for their valuable suggestions for this study.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Fig S1. Additional characteristics of Lund immunemethylation clusters. Fig S2. Additional characteristics of immune exclusion in the immune-methylation clusters. Fig S3. Additional characteristics of TCGA immunemethylation clusters. Table S1. List of immune cell type associated CpGs with corresponding gene and associated cell type information. Table S2. Clinical characteristics for immune-methylation clusters (Lund cohort). Table S3. Overlap between Lund gene expression phenotypes and Lund immune-methylation clusters. Table S4. List of centroid CpGs with methylation b values for each cluster. Table S5. Overlap between TCGA gene expression classes and TCGA immune-methylation clusters. Appendix S1. Supplementary methods.