DNA methylation–based immune response signature improves patient diagnosis in multiple cancers

Breast cancer (BC) remains challenging to treat because of its vast heterogeneous nature. Even within BC subtypes, patients experience different rates of survival and responses to anticancer therapies. This diversity needs to be further explored to improve prognosis and optimize therapeutic approaches in the future. In this respect, the tumor immune response is increasingly recognized to predict clinical outcomes in breast and other cancers (1). In particular, tumor-infiltrating lymphocytes (TILs) have been shown to provide prognostic and predictive information. High numbers of TILs have been associated with increased survival and response rates to preoperative chemotherapy in triple-negative (TN) (2–5) as well as human EGFR 2–positive (HER2) BCs treated with chemotherapy and trastuzumab (2, 4–6). Among TILs, T cell subsets are most abundant (7) and are associated with clinical outcomes (8–11) that suggest they play a major role in the antitumor immune response.

The evaluation of tumor immune responses by measuring TILs remains suboptimal, since pathologists base their quantification on subjective measurements. Microscopic counting of TILs using H&E- or IHC-stained tumor sections suffers from bias and variability and is only of a semiquantitative nature (12, 13). Recently, guidelines were published for a more consistent and reproducible morphological measurement of TILs, with the overall aim of establishing an “immunological grade” for clinical practice (14). Gene expression–based immune markers and signatures have been associated with TILs, and they also predict better clinical outcomes and response to therapy in TN and HER2 tumors (2, 5, 15).

DNA methylation plays a critical role in cell lineage specification and may therefore serve as a specific molecular marker for immune cell typing. Upon differentiation, cell lineage–specific changes occur in methylation, influencing the expression of key transcription factors and regulatory genes that lock the identity of cells (16–18). Although cell identity is determined by both DNA methylation and gene expression, DNA methylation may reflect distributions of cell subtypes more adequately, given that the relationship of only 2 DNA molecules per cell is of a more linear nature than are thousands of mRNA copies exposed to mRNA degradation (13, 19). Indeed, several studies have identified DNA methylation signatures that accurately evaluate the distribution of cell subpopulations in blood (19–21). However, within complex tissues such as tumors, DNA methylation has barely been explored in terms of evaluating immune cell subtypes, particularly TILs, and the diagnostic value of such a marker is unknown.

DNA methylation landscapes in BC are highly abnormal, and numerous studies have shown differences in BC methylomes according to clinical and pathological parameters such as the expression of hormone receptors, tumor grade and stage, and survival (22–25). Interestingly, differences captured through whole-tumor DNA methylome profiling not only originate from tumor cells but also from TILs, suggesting the potential of DNA methylation for the evaluation of tumor immune responses (26).

In this study, we applied genome-wide DNA methylation profiling to identify markers (methylation of TIL [MeTIL] signature) for evaluation of the local tumor immune response and their potential to improve the prognostic accuracy for BC patients. Our results showed that the MeTIL signature score measured TILs in a sensitive manner, resulting in a better prediction of survival for BC subtypes. We further showed that the MeTIL score might predict the response to anthracycline-based chemotherapy. This signature further stratifies patients with different clinical outcomes in other cancers including types for which TILs were not known to have a prognostic value. Finally, we demonstrated the clinical merit of applying this methodology in the clinic, since the MeTIL score can be determined by bisulfite pyrosequencing of low amounts of DNA from formalin-fixed, paraffin-embedded (FFPE) tumor tissue.

Evaluation of TIL distributions in breast tumors with the MeTIL signature. The MeTIL signature was developed in 2 steps and subsequently tested for its prognostic and predictive value in various cohorts, as outlined in Figure 1. Given that among TILs, T cells are the most abundant (7) and are associated with clinical outcomes (8–11), we identified, in a first step, CpGs that are highly differentially methylated between normal or cancerous breast epithelial cells and T lymphocytes (Supplemental Figure 1 and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI91095DS1). In a second step, we used DNA methylation profiles from breast primary tumors in cohort 1. Cohort 1 (n = 118) is an in-house retrospective cohort of BC patients who received adjuvant therapies according to institutional recommendations and were assessed by histopathological methods for the percentage of TILs, tumor cells, and other cell types of the tumor microenvironment (Supplemental Table 2). Pathological assessment of TILs (PaTILs) was performed on H&E-stained tumor sections by defining the percentage of mononuclear cells within the epithelium of the invasive tumor cell nests (Supplemental Table 3). We applied a random forest machine learning approach to DNA methylation profiles of cohort 1 to select markers from our list of 29 T cell–associated CpGs that most accurately predict the quantity of PaTILs in patients’ samples (Supplemental Figures 1 and 2). The final signature, named MeTIL, included 5 CpGs located within the promoter of 5 individual genes, namely protein tyrosine phosphatase, receptor type C–associated protein (PTPRCAP), internexin neuronal intermediate filament protein α (INA), semaphorin 3B (SEMA3B), Kelch-like family member 6 (KLHL6), and Ras association domain family member 1 (RASSF1) (Supplemental Table 4). As expected, their biological functions comprised immunity-related mechanisms (Supplemental Figure 1). The development of the MeTIL signature is described in detail in the Supplemental Methods.

Figure 1

Flow diagram of the cohorts used in the study. Cohort 1 and cohort 2 are in-house retrospective cohorts of patients diagnosed with luminal (LUM), HER2-positive, or TN breast tumors and treated with adjuvant therapies according to institutional recommendations. The TOP cohort consists of 58 samples from a prospective clinical trial (TOP), in which hormone receptor–negative BC patients were treated with preoperative epirubicin monotherapy. Note that 38 tumors in the TOP cohort were TN and used for the second validation of the prognostic value of the MeTIL score in TN tumors.

The ability of the MeTIL signature to evaluate TILs is based on highly differential methylation values of MeTIL markers in T cells, the major cell subtype among TILs, versus epithelial breast cells (Figure 2A). MeTIL markers also displayed differential methylation values in CD45⁺Epcam^– (lymphocytes) and CD45^–Epcam⁺ (epithelial) cells sorted from whole breast tumor samples (Figure 2B), suggesting that this signature discerns TILs from breast epithelial cells in vivo. Hierarchical clustering analysis based on the methylation values of MeTIL markers in cohort 1 identified 3 clusters with different methylation profiles (hypomethylated, intermediate methylated, and hypermethylated), which were associated with variable levels of PaTILs (P < 0.001) and distinct molecular subtypes (P = 0.003) (Figure 2C). This highlights the potential of the MeTIL signature to stratify tumors according to TILs.

Figure 2

Measurement of TIL distributions with DNA methylation (MeTIL signature). (A) Markers of the MeTIL signature showed highly differential methylation values in normal and cancerous breast epithelial cells (MCF10A, MCF-7, T47D, SKBR3, BT20, MDA-MB-231, MDA-MB-361, ZR-75-1) versus T lymphocytes (WEIS3E5, R12C9, and ex vivo T cells). (B) Unsupervised hierarchical clustering analysis of MeTIL marker β values in CD45⁺Epcam^– (lymphocytes) and CD45^–Epcam⁺ (epithelial) cells sorted from whole breast tumor samples. (C) Unsupervised hierarchical clustering analysis of cohort 1 based on β values of MeTIL markers. Note, a hypomethylated, an intermediate methylated, and a hypermethylated cluster appeared, all of which are associated with differences in BC subtypes, PaTILs, and MeTIL scores. Differences between methylation clusters were assessed with a 1-way ANOVA or χ² test, and P values are shown in the upper right corner of the heatmap. (D) The MeTIL signature was transformed into a score, and MeTIL scores were computed for T cells, Tregs, B cells, NK cells, granulocytes, monocytes, and DCs. Infinium DNA methylation profiles for these sorted blood cell populations are publicly available in the NCBI’s GEO database (GEO GSE35069, GSE39981, GSE49667, and GSE59796), per the Methods section. The MeTIL score values are plotted on the y axis and blood cell subpopulations on the x axis. Differences in the MeTIL score between the groups including T cells, Tregs, B cells, and NK cells and the group of granulocytes, monocytes, and DCs were assessed with a Student’s t test (*P < 0.05). (E) MeTIL scores were computed in tumors enriched for CTLs with high expression of GZMB and PRF1 GZMB/PRF1-high and in tumors enriched for CTLs with low expression of GZMB and PRF1 GZMB/PRF1-low. The difference in MeTIL scores between the 2 groups was assessed with a Student’s t test, and the P value is shown. (F) MeTIL scores were correlated with the percentage of adipocytes, fibroblasts, and endothelial cells or PaTILs for 62 samples from cohort 1. MeTIL score values are plotted on the y axis and the percentage of cells on the x axis. The correlation was assessed with a Spearman’s rank correlation test. The Spearman’s rank correlation coefficient (rho) and its P value are shown for each plot. Note, because of methodological limitations, fibroblasts and endothelial cells were assessed as 1 cell type. (G) Color map showing the MeTIL score performance for simulations across noise (y axis) and the presence of additional cell type(s) (x axis) (randomly selected methylation values). With a SD of noise of 1, the performance (assessed using the Spearman’s correlation metric R) stayed higher than 0.7, even if the tissue consisted of more than 70% non–BC cells and non–T cells (yellow border). (H) Cohort 1 (105 samples) and cohort 2 (100 samples) were grouped by PaTILs, and MeTIL scores are shown by PaTIL group. MeTIL score values are plotted on the y axis and PaTIL groups on the x axis. Differences in the MeTIL score between groups were assessed with a 1-way ANOVA, and the P value is shown in each plot. (I) TIL distributions in BC subtypes based on MeTIL score (left panel) or PaTILs (right panel) in cohort 1 and cohort 2. Note that the MeTIL score showed greater differences within (especially in LumA and LumB) and between subtypes than did PaTILs. BC subtypes were defined on the basis of IHC results for the hormone receptors and HER2.

To evaluate TIL distributions with the MeTIL signature, we transformed the individual methylation values of the MeTIL markers into a score (MeTIL score) by applying a principal component analysis (PCA) (described in the Supplemental Methods). We first sought to test which immune cell subtypes this score specifically measures. We used publicly available Infinium DNA methylation profiles from various sorted blood cell populations and computed the MeTIL score for T cells, Tregs, B cells, NK cells, granulocytes, monocytes, and DCs (Figure 2D) (27–30). Not surprisingly, T cells displayed the highest median MeTIL score. MeTIL scores were also high for Tregs, B cells, and NK cells, suggesting that these cell types contributed to the MeTIL score as well. Concordantly, MeTIL scores correlated strongly with total leucocytes (CD45⁺, rho = 0.506; P < 0.001), T cells (CD3⁺, rho = 0.651; P < 0.001), and B cells (CD20⁺, rho = 0.482; P < 0.001) in 61 samples from cohort 1 for which IHC staining data were available (Supplemental Figure 3). Granulocytes, monocytes, and DCs had significantly lower median MeTIL scores and are therefore less reflected by the MeTIL score (Figure 2D). Together, these results showed that the MeTIL score measures predominantly mononuclear immune cells including lymphocytes and NK cells and that among these cell subtypes, T cells contributed most to the MeTIL score. We further tested whether MeTIL scores displayed differences between T cell subsets. CD4⁺ and CD8⁺ cells sorted from blood showed no differences in the median MeTIL score (P = 0.130), suggesting that the MeTIL score reflects both helper and cytotoxic T lymphocytes (CTLs) (Supplemental Figure 4). Finally, we sought to test whether the MeTIL score reflects the functional states of CTLs. We used previously published expression data for cohort 1 (26) and defined tumors enriched for functional CTLs or nonfunctional CTLs on the basis of the expression of CD8A, granzyme B (GZMB) and perforin 1 (PRF1). CD8A is a marker for CTLs, and GZMB and PRF1 are both markers of CTL activity. Tumors with high expression of CD8A and high expression of both GZMB and PRF1 were therefore considered enriched for functional CTLs. Tumors with high expression of CD8A but low expression of both GZMB and PRF1 were considered enriched for nonfunctional CTLs (Supplemental Figure 5A). We found a significantly higher median MeTIL score in tumors enriched for functional CTLs (Figure 2E), suggesting that the MeTIL score is a measure of the functionality of CTLs. Notably, there were no differences in CD8A expression between tumors enriched for functional or nonfunctional CTLs, suggesting that the observed difference in MeTIL scores was not due to differences in CTL abundance but rather to differences in CTL functionality (Supplemental Figure 5B).

Next, we sought to test whether the MeTIL score is specific to TILs or whether it also reflects other cell types of the tumor microenvironment that have been shown to impact tumor progression and patient outcomes (31–33). Using Infinium DNA methylation data from the Encyclopedia of DNA Elements (ENCODE), we found the highest MeTIL scores for lymphocytes as compared with those for epithelial cells, fibroblasts, muscle cells, and other microenvironmental components such as adipocytes and endothelial cells (Supplemental Figure 6). We further correlated MeTIL scores with pathological assessments of adipocytes or fibroblasts and endothelial cells in 62 samples from cohort 1 and observed no significant correlation between MeTIL scores and adipocytes (rho = 0.083; P = 0.521) or between MeTIL scores and fibroblasts and endothelial cells (rho = –0.239; P = 0.061) (Figure 2F). In contrast, MeTIL scores strongly correlated with the PaTILs (rho = 0.696; P < 0.001) in these samples. Additionally, we estimated the MeTIL score performance in simulation models reflecting biological admixtures of breast tumors (34). The MeTIL score accurately resolved known mixture proportions with an unknown cell content of up to approximately 70% and a noise up to approximately 70% of the SD (Figure 2G). Together, these data suggest that the MeTIL score is specific to TILs and does not measure other cell types typically found in the breast tumor microenvironment.

After we confirmed that the MeTIL score specifically measures TILs, we tested this score for the evaluation of TIL distributions within breast tumors. When MeTIL scores were grouped by PaTILs, we observed a significant increase in median MeTIL scores, with increasing levels of PaTILs in cohort 1 (P < 0.001) and cohort 2 (P < 0.001) (Figure 2H), suggesting that MeTIL scores may be suitable for quantifying TILs. Cohort 2 (n = 119) is an in-house retrospective BC cohort that received adjuvant therapies according to institutional recommendations and was assessed for PaTILs (Supplemental Tables 3 and 5). We then applied the MeTIL score to measure TIL distributions according to BC subtypes, namely TN, HER2, luminal A (LumA), and luminal B (LumB). In cohorts 1 and 2, the MeTIL scores had the highest values in TN and HER2 tumors, corresponding to previous findings of the highest TIL abundance in these subtypes (Figure 2I). Interestingly, the MeTIL score showed wider distributions within and greater differences between subtypes than did PaTILs.

Improved prediction of survival and response to chemotherapy with the MeTIL score. Recent studies have shown that TILs carry prognostic information mainly in TN breast tumors (3, 4). We tested whether the MeTIL score and PaTILs predict survival within BC subtypes in cohort 1 and cohort 2. Additionally, we used 38 TN tumors from a previously published prospective clinical trial (Trial of Principle [TOP]), in which 149 hormone receptor–negative patients with BC were treated with preoperative epirubicin monotherapy (35) (Supplemental Table 6). Of note, PaTILs were not available for the TOP cohort. Clinical and pathological characteristics for all cohorts are shown in Supplemental Table 7. In TN tumors, high MeTIL scores, but not PaTILs, were shown to be associated with better outcomes in cohort 2 (hazard ratio [HR], 0.42; 95% CI, 0.20–0.91; P = 0.027) and in TN samples from the TOP trial (HR, 0.59; 95% CI, 0.39–0.89; P = 0.012) by univariate Cox analyses (Figure 3A and Supplemental Table 8). In luminal tumors, a high MeTIL score was associated with favorable outcomes in cohort 1 (HR, 0.65; 95% CI, 0.43–0.98; P = 0.041) and cohort 2 (HR, 0.25; 95% CI, 0.06–0.96; P = 0.043). These HRs remained borderline significant after correction for LumA/LumB status. HER2-positive tumors with high MeTIL scores showed a statistically nonsignificant trend toward better outcomes that was likely due to low sample numbers. Multivariate Cox analysis demonstrated that the prognostic value of the MeTIL score in TN tumors was independent of other clinical and pathological variables in the TOP cohort (HR, 0.43; 95% CI, 0.25–0.74; P = 0.003) (Table 1). In luminal tumors, we found borderline significant independent prognostic values for the MeTIL score in cohort 1 (HR, 0.69; 95% CI, 0.47–1.01; P = 0.055) and cohort 2 (HR, 0.17; 95% CI, 0.03–1.05; P = 0.056). Together, these data suggest that the MeTIL score may predict outcomes for TN tumors and tumors with low TILs infiltration (Lum).

Figure 3

The MeTIL score improves the prediction of survival and response to anthracycline treatment. (A) Forest plots by BC subtype showing the log₂ value of the HRs and CI for the prediction of survival outcomes in univariate Cox models for the MeTIL score (orange) or PaTILs (black) in 3 BC cohorts. BC subtypes were defined on the basis of IHC results for the hormone receptors and HER2. PaTILs were not available for the TOP cohort. Only ER- and HER2-negative tumors were selected from the TOP cohort as TN tumors. The red asterisk indicates statistical significance (P < 0.05 by a likelihood ratio test). (B) Receiver operating characteristic (ROC) curve for the prediction of response to neoadjuvant anthracycline treatment based on the MeTIL score for 58 hormone receptor–negative patients in the TOP cohort. (C) Forest plot showing the log₂ value for the OR and CI of the MeTIL score (orange) and various other clinical and pathologically relevant variables (black) for the prediction of response to preoperative anthracycline treatment in a multivariate analysis of the TOP cohort. The red asterisk indicates statistical significance (P < 0.05, which corresponds to the z ratio based on a normal reference distribution).

Table 1

Correlation between the MeTIL score or PaTILs and median survival of BC patients in the context of other prognostic clinical and pathological variables by BC subtype in various cohorts (multivariate Cox proportional hazards regression)

TILs have been associated with higher response rates to preoperative chemotherapy in hormone receptor–negative and HER2-positive BCs (5). We assessed the potential of the MeTIL score to predict the response to preoperative chemotherapy in 58 hormone receptor–negative breast tumors from the TOP trial. An AUC of 0.73 (95% CI, 0.547–0.92) suggested a predictive value for the MeTIL score (Figure 3B). Logistic regression modeling demonstrated that the MeTIL score predicts for response to preoperative chemotherapy independently of other clinical and pathological variables, with an odds ratio (OR) of 4.38 (CI, 1.62–21.5; P = 0.02) (Figure 3C and Supplemental Table 9). Of interest, the MeTIL score showed the highest AUC and a significant OR (Supplemental Figure 7, A and B, and Supplemental Table 10) for response when compared with several gene expression–based immune markers (CD3D, CXCL9, CD247) (15) and signatures such as the STAT1 metagene (36). These results suggest, if further validated, that the MeTIL score may be a potential marker of response to chemotherapy in the future.

Evaluation of TILs in low amounts of DNA from FFPE tumor tissue through bisulfite pyrosequencing of MeTIL markers. The MeTIL score measures TILs in a sensitive manner and results in an improved prediction of survival and response to therapy. It would therefore be an attractive tool in clinical practice. To make its application easy, fast, and feasible in the clinic, we optimized bisulfite pyrosequencing for MeTIL score measurement in FFPE tumor tissue. We sequenced 21 FFPE tumor samples from cohort 1, for which Infinium array–based MeTIL scores from frozen tissue were available, and observed significant correlations between the methylation values of MeTIL markers obtained by bisulfite pyrosequencing (y axis, Figure 4A) and Infinium arrays (x axis, Figure 4A). MeTIL scores obtained through bisulfite pyrosequencing of individual markers strongly correlated (rho = 0.79, P < 0.01) with Infinium-based MeTIL scores (Figure 4B). Of note, bisulfite pyrosequencing–based MeTIL scores stratified breast tumors according to PaTILs (Figure 4B) and subtypes, since these are associated with different levels of TILs (Figure 4C).

Figure 4

Determination of the MeTIL score in FFPE tumor tissue by bisulfite pyrosequencing. (A) Scatter plots showing for each MeTIL marker the correlation between methylation values (percentage), determined by pyrosequencing in FFPE tumor tissue (y axis), and methylation values (percentage), determined by Infinium arrays in fresh-frozen tumor tissue (x axis). The correlation was established on the basis of 21 samples from cohort 1, for which fresh-frozen and FFPE tissue was available. Spearman’s rank correlation coefficient (rho) and its P value for each marker are shown. Different colors of the dots reflect the PaTIL group to which each sample was assigned. PaTIL groups were defined on the basis of TIL percentages as follows: no PaTILs (PaTILs <1%); low PaTILs (PaTILs ≥1% and <20%); and high PaTILs (PaTILs ≥21% and ≤100%). (B and C) Scatter plots showing the correlation between the MeTIL score, determined by pyrosequencing in FFPE tumor tissue (y axis), and the MeTIL score, determined by Infinium arrays in fresh-frozen tumor tissue (x axis), with respect to the PaTIL group (B) or BC subtype (C).

Prediction of survival outcomes in other cancer types with the MeTIL score. Acknowledging that the MeTIL signature was developed for the evaluation of TILs in BC, we assessed whether it could predict survival differences in other cancer types available in TCGA (Supplemental Table 11). In 5 of 20 tested cancer types (head and neck squamous cell carcinoma [HNSC], pheochromocytoma and paraganglioma [PCPG], skin cutaneous melanoma [SKCM], thyroid carcinoma [THCA], and thymoma [THYM]), high MeTIL scores, but not PaTILs, were associated with a better outcome (Figure 5A and Supplemental Table 12). MeTIL scores predicted survival differences independently of other prognostic variables in HNSC, PCPG, SKCM, and THYM, but also in lung squamous cell carcinoma (LUSC) (Figure 5B and Supplemental Table 13). As with BC (Figure 2C), MeTIL scores clustered SKCM samples into 3 groups that were associated with variable levels of PaTILs (P = 0.049), distinct molecular subtypes (37) (P < 0.001), and MeTIL scores (P < 0.001) (Figure 5C) as well as differences in survival (P = 0.018) (Figure 5D). MeTIL scores varied for melanoma subtypes, with the “immune” subtype showing the highest median MeTIL score (P < 0.001) (Figure 5E). Together, these data suggest that the MeTIL score may predict survival outcomes in other cancer types. It is noteworthy that we did not observe significant differences in BC survival with the MeTIL score in TCGA data. This was not unexpected, as TILs have been shown to be more abundant and associated with clinical outcomes in TN and HER2 breast tumors. We therefore correlated MeTIL scores with survival endpoints in BC subtypes and observed differences in survival with HER2 tumors (HR, 0.37; 95% CI, 0.16–0.85; P = 0.02). In luminal (HR, 0.82; 95% CI, 0.67–1.02; P = 0.069) and TN (HR, 0.67; 95% CI, 0.44–1.03; P = 0.066) tumors, the association between MeTIL score and survival was borderline significant. PaTILs, on the other hand, predicted no survival differences between luminal (HR, 0.98; 95% CI, 0.93–1.02; P = 0.324), HER2 (HR, 0.56; 95% CI, 0.20–1.58; P = 0.273), or TN (HR, 0.99; 95% CI, 0.96–1.03; P = 0.733) tumors. These results are in line with our earlier findings that suggested a prognostic value for the MeTIL score in luminal and TN tumors and further show a prognostic value for HER2 tumors. Importantly, as in BC, the MeTIL score may have a prognostic value in other cancer types if these are grouped into subtypes.

Figure 5

The MeTIL score predicts differences in survival in other types of cancer. Forest plot showing the log₂ value of the HR and CI for the prediction of survival outcomes in univariate (A) or multivariate (B) Cox models for the MeTIL score (orange) or PaTILs (black) in different TCGA cancer types. Red asterisks indicate statistical significance (P < 0.05 by a likelihood ratio test). BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; COREAD, colon and rectum adenocarcinoma; ESCA, esophageal carcinoma; KIRP, kidney renal papillary cell carcinoma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PRAD, prostate adenocarcinoma; SARC, sarcoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; UCEC, uterine corpus endometrial carcinoma. (C) Heatmap displaying the results of an unsupervised hierarchical clustering analysis of TCGA skin cutaneous melanomas based on β values for the MeTIL markers. Note, a hypomethylated, an intermediate methylated, and a hypermethylated cluster appeared, all of which are associated with differences in subtypes, PaTILs, and MeTIL scores. Differences between methylation clusters were assessed by 1-way ANOVA (MeTILs) or χ² test (PaTILs and subtypes), and P values are shown. (D) Kaplan-Meier survival curves for the 3 methylation clusters defined in the heatmap. (E) MeTIL scores grouped according to 3 melanoma subtypes. Differences in MeTIL scores between melanoma subtypes were assessed by 1-way ANOVA, and the P value is shown. MITF, malenogenesis-associated transcription factor.

The tumor immune response is increasingly recognized to be associated with better clinical outcomes in breast and other cancers. However, quantitative evaluation of the tumor immune response based on TILs remains suboptimal, since histopathological measurements are semiquantitative and limited in their accuracy and reproducibility (12, 13). This has prompted an international TIL working group to publish guidelines for the harmonization of this method (14). DNA methylation plays a critical role in cell lineage specification and may therefore sustain a specific molecular marker for typing of immune cell subtypes (16–18). Indeed, several studies have identified DNA methylation signatures that accurately evaluate the distribution of cell subpopulations in blood (19–21). However, within complex tissues such as tumors, DNA methylation has barely been explored for the evaluation of immune cell subtypes, particularly TILs, and the diagnostic value of such a marker is unknown.

In this study, we identified DNA methylation markers (MeTIL signature) for the evaluation of TIL-based tumor immune responses and their impact on clinical outcomes in BC. Interestingly, the functions of all MeTIL markers were related to TILs or other components of the antitumor immune response. PTPRCAP is a phosphoprotein that is specifically associated with CD45, a surface marker on lymphocytes, and has been shown to function as a key regulator of T and B cell activation (38). KLHL6 is a member of the Kelch-like protein family and is important for antigen receptor signaling on B cells and germinal center formation (39). INA and SEMA3B are best known for their role in neuronal development, but studies also linked them to immune-related functions (40, 41). INA encodes for the neurofilament protein internexin neuronal intermediate filament protein α, and neurofilaments have been linked to T cell activation (40). SEMA3B belongs to the family of semaphorins, which regulate immune functions by controlling activation, differentiation, and trafficking of immune cells, including T cells and B cells (41). RASSF1 is a well-established tumor suppressor that is frequently inactivated in several cancers, including breast and lung cancer, by aberrant promoter methylation. RASSF1 controls genome stability in response to replication stress through activation of the Hippo pathway, which regulates phosphorylation of breast cancer 2 (BRCA2) and recruitment of RAD51 recombinase (RAD51) (42). The absence of RASSF1A led to chromosomal aberrations and increased genomic instability such as that seen in BRCA-mutant cells (42). Concordantly, a positive correlation between RASSF1A promoter methylation and increased copy number alteration has been shown in breast and lung cancer (43). These findings, together with those of other studies, which linked genomic instability to more TILs and a better antitumor immune response, suggest that methylation of RASSF1 may be an indirect measure of TILs and the antitumor immune response (44, 45).

A thorough characterization of the MeTIL signature showed that the signature score measures predominantly mononuclear immune cells including T cells, B cells, and NK cells. Among these cell subtypes, T cells had the highest MeTIL scores, suggesting a bias of the MeTIL score toward T cells. Moreover, MeTIL scores were markedly higher in tumors enriched for functional CTLs than in those enriched for nonfunctional CTLs, suggesting that the MeTIL score may reflect the functionality of tumor immune responses. Intriguingly, MeTIL scores were also high in Tregs, which have an immunosuppressive role in tumors. This is interesting in light of previous studies that showed a positive correlation between immunosuppressive markers, including the Treg marker FOXP3, and TILs and suggested a feedback activation of immunosuppressive pathways as part of the immune reaction (5, 46). Hence, Tregs and other immunosuppressive markers can be a surrogate for an immune reaction in tumors, and their contribution to the MeTIL score further facilitates the quantification of tumor immune responses. Other immune cells, e.g., granulocytes, monocytes, and DCs, showed markedly lower MeTIL scores and hence contributed only minorly to the MeTIL score. Importantly, the frequencies of nonimmune cells (adipocytes, fibroblasts, endothelial cells), which are typically found in the tumor microenvironment, did not correlate with MeTIL scores. Together, these results suggest that the MeTIL score measures predominantly TILs. This is substantiated by the strong correlation we observed between the MeTIL score and PaTILs in vivo. Next, we used the MeTIL score to measure TIL distributions within breast tumors of different subtypes in 2 independent cohorts. We consistently observed differences between subtypes with the MeTIL score, but not with PaTILs. In line with other studies, the MeTIL score showed the highest TIL levels in TN and HER2 tumors (3, 4). Interestingly, also within BC subtypes, especially in luminal tumors with low infiltration, the MeTIL score showed wider TIL distributions than did PaTILs. Together, these results may suggest a greater sensitivity of the MeTIL score for the evaluation of TILs. The higher sensitivity of the MeTIL score may be attributed to the specificity and stability of DNA methylation markers for TIL detection as well as the sensitive array technology used to measure these marks. The greater tumor area used for methylation analysis than for pathological assessment (typically 1 slide per tumor sample) (14) may further explain the discrepancy between MeTIL score–based and pathology-based TIL levels.

TILs have recently been reported to predict survival outcomes, particularly in TN and HER2-positive BCs (3, 4, 47). We therefore evaluated the MeTIL score for its accuracy in predicting survival within BC subtypes and found an association between high MeTIL scores, but not high PaTILs, and better clinical outcomes for TN and luminal tumors. To our knowledge, this is the first time a prognostic value of immune-related markers has been demonstrated in luminal tumors. Luminal tumors are generally less infiltrated by lymphocytes, making it challenging to establish an association with outcome. However, the high sensitivity of the MeTIL score for the TIL-based tumor immune response and, potentially, its bias toward T cells and functional CTLs may have allowed us to observe such an association. We further showed that high MeTIL scores were associated with higher response rates to preoperative chemotherapy in patients with hormone receptor–negative BC. The MeTIL score remained the only significant marker in a multivariate analysis and performed at least as well as expression-based immune markers. Interestingly, in our analyses, PaTILs were not associated with clinical outcome in TN tumors, although other studies have established PaTILs as a prognostic marker in this subtype (3, 4). This discordance may be explained by the much smaller cohort sizes we used in our study as compared with the large sample numbers that were used in other studies to reveal the prognostic value of PaTILs. Nevertheless, the cohort sizes in our study were sufficiently powerful to reveal an association between MeTIL scores and survival or response to therapy.

Given the sensitive nature of the MeTIL score for TIL evaluation as well as its improved prediction of survival outcome and response to therapy as compared with, respectively, PaTILs or gene expression–based immune markers, the MeTIL score could be an attractive tool for evaluating immune responses in clinical settings. To make this feasible, we optimized bisulfite pyrosequencing as a readout of the MeTIL score from FFPE tumor tissue and showed that this MeTIL score strongly correlated with the MeTIL score generated from Infinium methylation arrays from frozen tumor tissue. Thus, bisulfite pyrosequencing of MeTIL markers may be used as a fast, easy, and cost-efficient method for TIL evaluation in the clinic.

Intriguingly, the MeTIL score was able to improve the accuracy of prognosis in other types of cancer. In SKCM, LUSC, HNSC, PCPG, THCA, and THYM, high MeTIL scores, but not high PaTILs, were associated with better outcomes. While the prognostic value of PaTILs has been well established in melanoma (48), lung cancer (49), and HNSC (50), the MeTIL score highlights, for the first time to our knowledge, a prognostic value related to TILs in PCPG, THCA, and THYM. Together, these results suggest that the MeTIL signature, although developed for the evaluation of tumor immune responses in breast tumors, may also measure immune responses in other types of cancer and thus more accurately stratify patients in terms of prognosis. Interestingly, the MeTIL score did not predict for survival differences in BC from TCGA data, but when tumors were grouped by subtype, a prognostic value was observed in HER2 tumors and, albeit to a lesser degree, in luminal and TN tumors. While these results support our earlier finding that revealed a prognostic value of the MeTIL score in luminal and TN tumors, they also suggest that the MeTIL score may have prognostic value in other cancer types from TCGA when grouped by subtype. Nevertheless, we were intrigued to find no association between PaTILs and clinical outcome in any of the cancer types, although PaTILs has been established as a prognostic marker in several of them (48–50). As discussed by Fridman et al., the prognostic value of PaTILs has been established in large clinical studies involving more than 1,000 patients each. Most TCGA cohorts in our study included substantially fewer patients, and this may have hampered a possible association between PaTILs and survival. Importantly, TILs have to be appropriately analyzed to provide prognostic information (12, 14). PaTILs suffers from bias and variability and is only semiquantitative in nature (12). Consequently, variability among PaTILs is generally high, and this may have masked the association between PaTILs and clinical outcomes. The MeTIL score on the other hand, measures TILs in an unbiased manner and thus is less prone to variability in TCGA cohorts. The unbiased evaluation of the MeTIL score allowed us to establish strong associations with survival, despite the smaller cohort sizes. The discrepancy in prognostic performance between the MeTIL score and PaTILs may be further explained by differences in the characteristic of the variable. The Immuno-Oncology Biomarker Working Group (previously known as the TIL Working Group) suggested that TILs should be scored as a continuous variable, as this allows for more accurate statistical analysis and thus more accurate associations with clinical variables such as prognosis (14). In daily practice, however, TILs are considered mostly in categories of 5% each, since the eye of the pathologist cannot distinguish between, for example, 11% of TILs and 14% of TILs and thus will score this example as either 10% or 15% (14). This categorical assessment is suboptimal, as relevant information, which is important for the prognostic TIL effect that has been shown to be linear in various studies, is lost (51). The guidelines therefore suggest that TILs be scored as a continuous variable and as accurately as possible, which means that TIL categories should be kept as small as possible to avoid the loss of prognostic and potentially predictive information. Nevertheless, in contrast to the MeTIL score, which is a real continuous parameter, PaTILs is, strictly speaking, rarely scored as a continuous variable in daily clinical practice (14). This probably caused some loss of relevant prognostic information with regard to PaTILs and might further explain discrepancies with MeTIL scores and with studies that have assessed TILs as a continuous variable. Last, we demonstrated that the MeTIL score predominantly reflects T cells and functional CTLs. This bias further adds to the discrepancy in prognostic performance between MeTIL scores and PaTILs.

In this study, we highlighted the power of DNA methylation to evaluate local and functional TIL-based tumor immune responses and the ability of this approach to improve prognosis in breast and other cancers. The MeTIL signature, if further validated, holds potential for the future to refine the stratification of cancer patients for clinical trials and the choice of therapeutic approaches, including immunotherapy.

Patient cohorts. Cohort 1 and cohort 2 consists of 118 and 119 retrospectively selected fresh-frozen tumor samples from patients treated with adjuvant therapies according to institutional recommendations and diagnosed at the Jules Bordet Institute from 1995 to 2003 and 2004 to 2009, respectively. The preoperative TOP cohort consisted of 149 patients with estrogen receptor–negative (ER-negative) disease who were treated at the Jules Bordet Institute from 2003 to 2008 with neoadjuvant epirubicin monotherapy (100 mg/m²) (35). Patients with operable BC were treated every 3 weeks for 4 cycles, and patients with inflammatory or locally advanced BC were treated every 2 weeks for 6 cycles. Pretreatment biopsies were obtained from the primary lesion. A pathologic complete response (pCR) was the primary endpoint of this trial. A pCR was defined as the absence of residual invasive breast carcinoma in the breast and in the axillary nodes after completion of chemotherapy. Persistence of in situ carcinoma without an invasive component was also considered a pCR. Fifty-eight fresh-frozen samples with a sufficient number of tumor cells and high DNA yield were selected for Infinium methylation analyses. The clinical and pathological characteristics of the 3 patient cohorts are summarized in Supplemental Table 7.

Histopathological assessment of TILs and other cell types of the tumor microenvironment. Histopathological assessment of TILs (PaTILs) was performed on H&E-stained tumor sections by defining the percentage of mononuclear cells within the epithelium of the invasive tumor cell nests (Supplemental Table 3). TIL frequencies were evaluated independently by 2 well-trained pathologists, and the mean value was used for analyses.

PaTIL results were available for cohort 1, cohort 2, and TCGA tumors, but not for the TOP cohort. In 61 samples from cohort 1, total leukocytes (CD45⁺), T cells (CD3⁺), and B cells (CD20⁺) were assessed by IHC (Supplemental Table 3). Assessments of other cell types in the tumor microenvironment were performed on H&E-stained tumor sections. The total percentage of lymphocytes, tumor cells, adipocytes, fibroblasts, and endothelial cells was scored on 1 full section per case. Because of methodological limitations, fibroblasts and endothelial cells were scored as 1 cell type.

TCGA and ENCODE data. DNA methylation and clinical data for various TCGA cancer cohorts were downloaded as described in the supplemental material. The clinical and pathological characteristics for TCGA cohorts are summarized in Supplemental Table 11. Histopathological measurements of TILs (PaTILs) for TCGA tumors were performed as summarized in Supplemental Table 3. ENCODE data for different cell types typically found in a breast tumor biopsy (breast epithelial cells, lymphocytes, fibroblasts, muscle cells, endothelial cells, and adipocytes) were retrieved and processed as described in the supplemental material.

Tumor preparation and flow cell sorting. Fresh breast tumor tissues were collected immediately following surgery and dissociated (without enzymatic digestion) using the gentleMACS Dissociator (Miltenyi Biotec) prior to Ab labeling, as previously described (52). Fluorescence-conjugated Abs against CD45 (BD Biosciences) and Epcam (Miltenyi Biotec) were used for surface staining of cells, according to the manufacturer’s protocol. Fluorescently labeled lymphocytes (CD45⁺Epcam^–) and tumor cells (CD45^–Epcam⁺) were sorted on a MoFlo Astros EQ 12/4 cell sorter. Cell purity was controlled on a Gallios 10/3 Cytometer and analyzed using Kaluza 1.3 Flow Analysis Software (both from Beckman Coulter).

Infinium HumanMethylation450K. Genomic DNA was extracted with the DNeasy Blood and Tissue Kit or the QIAamp DNA Mini Kit (both from QIAGEN) as previously described (26). DNA methylation was analyzed on Infinium HumanMethylation450K bead arrays as previously described for Infinium HumanMethylation27K bead arrays (26). Briefly, genomic DNA (300–800 ng) was converted with sodium bisulfite using the EZ DNA Methylation Kit (Zymo Research), and methylation assays were performed with 4 ml converted DNA at 50 ng/ml, according to the manufacturer’s protocol. Infinium HumanMethylation450K raw data were submitted to the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE72308; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72308).

Bisulfite pyrosequencing. Genomic DNA (275 ng) was bisulfite converted with the EZ DNA Methylation Kit, and 3–6 μl converted DNA (corresponding to approximately 45 to 95 ng DNA) served as the input for PCR. PCR assays were performed with HotStarTaq DNA Polymerase (QIAGEN) under the following cycle conditions: 95°C for 15 minutes, 50–60 cycles at 95°C for 1 minute, 50 cycles at 53°C for 1 minute, 72°C for 1 minute, and 72°C for 10 minutes. Amplification was confirmed on agarose gel, and pyrosequencing of successfully amplified PCR products was performed with the PyroMark Q24 System (QIAGEN). Primer sequences are listed in Supplemental Table 14.

Bioinformatics. Infinium HumanMethylation450K raw data (uncorrected probe intensity values) were preprocessed and β values computed and corrected as described in the supplemental material. T lymphocyte–associated markers (Supplemental Table 1) were identified through an approach using previously published DNA methylation profiles from normal or cancerous breast epithelial cell lines and T lymphocyte samples (26) (Supplemental Figure 1). From these, MeTIL markers (Supplemental Table 4) were selected by applying machine learning to cohort 1 (Supplemental Figures 1 and 2). Individual methylation values of MeTIL signature markers were transformed into a score (MeTIL score) using a normalized PCA (NPCA) approach (Supplemental Table 15). To estimate the performance of the MeTIL score, various biological admixtures of solid tumors were simulated as reported by Newman et al. (34). A detailed description of the bioinformatic methods is provided in the supplemental material.

Statistics. Statistical analyses were conducted with RStudio, version 0.94.110. Differences between more than 2 groups were assessed with a 1-way ANOVA or χ² test. Cox proportional hazard regression analyses and Kaplan-Meier survival curves with log-rank tests, recording patients at the time of dead or disease recurrence or last follow-up visit, were used to compare overall survival or disease-free survival rates. Multivariate Cox regression models were established on the basis of Akaike’s information criterion (AIC). ORs were used to compare pCR rates. The AUC was used to assess prediction performance. All P values were 2 sided, and P values of less than 0.05 were considered statistically significant. Statistical methods are further explained in the Supplemental Methods.

Study approval. This study was approved by the Medical Ethics Committee of Institute Jules Bordet, Brussels, Belgium, and all patients gave written informed consent before their participation in the study.