ARTICLE https://doi.org/10.1038/s41467-019-12159-9 OPEN DNA methylation loss promotes immune evasion of tumours with high mutation and copy number load Hyunchul Jung 1,10, Hong Sook Kim2,10, Jeong Yeon Kim1, Jong-Mu Sun2, Jin Seok Ahn2, Myung-Ju Ahn2, Keunchil Park2, Manel Esteller3,4,5,6,7, Se-Hoon Lee2,8 & Jung Kyoon Choi1,9 1234567890():,; Mitotic cell division increases tumour mutation burden and copy number load, predictive markers of the clinical benefit of immunotherapy. Cell division correlates also with genomic demethylation involving methylation loss in late-replicating partial methylation domains. Here we find that immunomodulatory pathway genes are concentrated in these domains and transcriptionally repressed in demethylated tumours with CpG island promoter hyper- methylation. Global methylation loss correlated with immune evasion signatures indepen- dently of mutation burden and aneuploidy. Methylome data of our cohort (n = 60) and a published cohort (n = 81) in lung cancer and a melanoma cohort (n = 40) consistently demonstrated that genomic methylation alterations counteract the contribution of high mutation burden and increase immunotherapeutic resistance. Higher predictive power was observed for methylation loss than mutation burden. We also found that genomic hypo- methylation correlates with the immune escape signatures of aneuploid tumours. Hence, DNA methylation alterations implicate epigenetic modulation in precision immunotherapy. 1 Department of Bio and Brain Engineering, KAIST, Daejeon 34141, Republic of Korea. 2 Division of Hematology/Oncology, Department of Medicine, Samsung Medical Center, Sungkyunkwan University School of Medicine, Seoul 06351, Republic of Korea. 3 Cancer Epigenetics and Biology Program (PEBC), Bellvitge Biomedical Research Institute (IDIBELL), L’Hospitalet, Barcelona, Catalonia, Spain. 4 Centro de Investigacion Biomedica en Red Cancer (CIBERONC), 28029 Madrid, Spain. 5 Institucio Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Catalonia, Spain. 6 Physiological Sciences Department, School of Medicine and Health Sciences, University of Barcelona (UB), Barcelona, Catalonia, Spain. 7 Josep Carreras Leukaemia Research Institute (IJC), Badalona, Barcelona, Catalonia, Spain. 8 Department of Health Sciences and Technology, Samsung Advanced Institute of Health Science and Technology, Sungkyunkwan University, Seoul 06351, Republic of Korea. 9 Penta Medix Co., Ltd., Seongnam-si, Gyeongi-do 13449, Republic of Korea. 10These authors contributed equally: Hyunchul Jung, Hong Sook Kim. Correspondence and requests for materials should be addressed to S.-H.L. (email: shlee119@skku.edu) or to J.K.C. (email: jungkyoon@kaist.ac.kr) NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 C ancer immunotherapy based on checkpoint blockade has cancer type (Supplementary Fig. 1). Our measure of genomic become highly effective in a subset of patients with dif- demethylation based on long interspersed nuclear element-1 ferent types of human cancers. In particular, antibody- (LINE-1 or L1)24,25 probes (Supplementary Fig. 2) also strongly mediated interventions targeting cytotoxic T lymphocyte antigen- correlated with cell proliferation markers (Fig. 1a and Supple- 4 (CTLA-4) and programmed death receptor-1 (PD-1) on T mentary Data 1). Global methylation loss was also associated with lymphocytes and the principal ligand (PD-L1) on tumour cells an increase in mutation burden (Fig. 1b) and chromosomal can reverse tumour-induced immunosuppression and induce SCNA load (Fig. 1c), two types of genomic aberrations that durable clinical responses1. accumulate through cell division (Supplementary Fig. 5). A major challenge facing current immunotherapies is the Notably, we found a correlation between the global L1 identification of biomarkers that predict clinical responses to methylation levels and immune signatures such as markers of CTLA-4 and PD-1/PD-L1 blockade. Overall, the mutational or tumour-infiltrating CD8 + T cells (Fig. 1d, Supplementary Data 1, neoantigen load2–5 and pre-existing T cell infiltration6,7 are and Supplementary Fig. 3). However, immune cell markers are indicators of clinical benefit of checkpoint blockade. On the other expected to correlate with mutation burden and also are known to hand, somatic copy number alterations (SCNAs)8–10, tumour be associated with aneuploidy10. To disentangle this intercorrela- heterogeneity11, and the genetic alteration of specific genes12 or tion, we performed multiple regression of the expression level of pathways13 have been identified as resistance factors. each gene on sample-level features, namely, global L1 methyla- Tumour cells produce neoantigens or antigens that the tion, mutation burden, aneuploidy, tumour purity, age, and immune system never encountered without cancer. The epitopes tumour stage. In this manner, we were able to determine that of neoantigens are displayed on the surface of cancer cells and immune infiltrates are associated with the global methylation provoke immune response. Therefore, tumours with high muta- levels independently of mutation burden and aneuploidy when tion load are more likely to respond to anti-immunosuppressive purity, age, and tumour stage are adjusted (Fig. 2a). strategies based on checkpoint blockade2–5. Mutation load Significant correlations with genomic demethylation were increases as a result of replicative errors during cell division. Not observed also for immunomodulatory pathways that should only mutations but also methylation losses accumulate during include genes expressed in tumour cells. These include antigen successive rounds of cell division14. Global hypomethylation and processing and presentation, major histocompatibility complex CGI hypermethylation represent the hallmark methylation (MHC), cytokine–cytokine receptor interaction, interferon or changes in cancer15. Considering its association with late repli- other cytokine signaling, and complement and coagulation cation timing, progressive methylation loss may occur due to the (Fig. 2b). There is an emerging role for the complement system failure of methylation maintenance machinery to remethylate in regulating the antitumour immune response26. The correlation newly synthesized daughter stands during DNA replication14,16. of the cell proliferation markers was in the opposite direction to However, the influence of genomic methylation loss through cell that of the immune cell markers or immunomodulatory genes divisions has never been investigated in the context of cancer (Fig. 2b). We confirmed that the global L1 methylation level itself immunotherapy. was not affected by the leukocyte fraction, which correlated only Meanwhile, SCNAs emerged as a resistance parameter8–10. A with immune gene expression (Supplementary Fig. 4). pan-cancer analysis discovered the association of SCNAs with molecular signatures of cytotoxic immune activity across diverse tumour types10. Particularly, highly aneuploid tumours with Repression of immune genes in late-replicating regions. To extensive chromosome- or arm-level SCNAs showed a lower focus on the repression of immune-response genes in tumour expression of markers indicating infiltrating immune cells. In cells, we excluded genes that are specifically expressed in the contrast, focal SCNAs mainly correlated with cell proliferation immune system from the following analyses. Because methylation markers instead of immune activity signatures. However, the loss occurs primarily in late-replicating regions14,16, we examined mechanism by which aneuploidy affects immune cell infiltration whether the transcriptional activity of late-replicated genes are remains unknown. Global demethylation in cancer promotes affected in the tumours that underwent global demethylation. By chromosomal instability17–20, particularly involving large-scale using cell line data, we identified genes that are replicated earlier alterations leading to aneuploidy21–23. Therefore, we investigated or later in cancer compared with normal cells (Supplementary the relationships between methylation changes and aneuploidy. Data 2). As a result, we found that the genes replicating late in Here, we performed large-scale systematic analyses of the cancer were significantly repressed in the demethylated tumours molecular data of TCGA samples across a variety of tumour (Fig. 3a) with CpG island (CGI) promoter hypermethylation types. We examined the relationships of global methylation levels (Fig. 3b). In contrast, the early-replicating genes tended to be with markers of cell proliferation, mutation burden, SCNA levels, overexpressed in the demethylated tumours (Fig. 3c). markers of infiltrating immune cells, and activity of immune- Overall, immune-related pathways were overrepresented in the response genes. Importantly, we tested our hypotheses developed late-replicating regions while cell cycle genes were concentrated in from the molecular analyses by using our lung cancer cohort. the early-replicating regions (Fig. 3d). More specifically, the This is the first study that inspected DNA methylation patterns in pathways most enriched for the late-replicating genes in cancer the molecular and clinical data with regards to cancer immu- included cytokine-cytokine receptor interaction, interferon-α/β notherapy. As a result, we suggest that as an important predictive (IFN-α/β) signaling, and RIG-1/MDA5-mediated IFN-α/β induc- marker in immunotherapy, genomic demethylation implicates tion (Fig. 3e and Supplementary Table 1). RIG-1/MDA5-mediated epigenetic modulation as a combination regimen for precision induction of IFN-α/β represents innate immune reaction against immunotherapy. RNA viruses. In contrast to our data from tumour-intrinsic demethylation, treatment of methylation inhibitors was shown to induce double-stranded RNAs (dsRNAs) derived from endogen- Results ous retroviruses (ERVs) and LINEs, resulting in the activation of Global methylation correlates with immune signatures. Our the IFN-α/β response in cancer27–29. Without the silencing of the pan-cancer analyses of TCGA data demonstrated that markers of IFN-α/β pathway, genomic demethylation would cause the cell proliferation tightly correlate with mutation burden and antiviral response and facilitate antitumour immune reaction as aneuploidy across cancer types and among samples within each demonstrated with demethylating agents. We measured the 2 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 ARTICLE a b r = –0.57 CESC r = –0.51 1000 TGCT P = 1.2 × 10 –2 P = 6.7 × 10–3 SKCM LUSC HNSC 0.2 BLCA CRC LUSC SKCM BLCA ESCA SARC LUAD Proliferation score ESCA Mutation burden STAD CRC STAD BRCA HNSC CESC LIHC 0.1 100 UCEC KIRP LUAD KIRC UCEC SARC LIHC BRCA PAAD PAAD PRAD LGG 0.0 TGCT KIRC PRAD LGG THCA KIRP PCPG THCA PCPG 10 0.35 0.40 0.45 0.50 0.55 0.35 0.40 0.45 0.50 0.55 Global methylation level Global methylation level Partial correlation coefficient (PCC) 0.0 0.0 –0.2 PCC –0.2 –0.4 –0.4 BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC c d r = 0.62 TGCT r = –0.56 KIRC –3 –3 P = 2.7 × 10 P = 8.6 × 10 0.22 LUSC PAAD 15 ESCA LUAD BLCA PRAD LGG Aneuploidy level CRC CD8 T cells SKCM HNSC 0.21 THCA LUAD 10 BRCA SARC STAD SARC BRCA LIHC CESC KIRP PAAD PCPG TGCT 5 BLCA 0.20 CRC PCPG KIRC KIRP UCEC SKCM CESC LGG LUSC PRAD LIHC 0 THCA HNSC 0.35 0.40 0.45 0.50 0.55 0.35 0.40 0.45 0.50 0.55 Global methylation level Global methylation level 0.0 0.4 PCC PCC –0.2 0.2 –0.4 0.0 BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC Fig. 1 Correlates of global DNA methylation level. Correlation between genomic methylation levels and cell proliferation markers a, mutation burden b, aneuploidy level c, and tumour-infiltrating CD8 + T cell markers d across and within 21 cancer types. a–d The median values were obtained for each cancer type, and statistical significance was evaluated using Spearman’s correlation (upper scatterplots). Three outlier cancer types (ESCA, STAD, and UCEC) are not shown for the CD8 + T cell correlation, but included when evaluating Spearman’s correlation. For the correlation within each cancer type (lower bar graphs), Spearman’s partial correlation was used to adjust for tumour purity. Tumour types showing significant partial correlation coefficient (P < 0.05) were shaded in darker colours NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 a Predictor variables (purity, age, and stage adjusted) Global methylation level Mutation burden Aneuploidy level Tumor type Activated_B_cell Activated_CD4_T_cell Activated_CD8_T_cell Central_memory_CD4_T_cell Central_memory_CD8_T_cell Effector_memeory_CD4_T_cell Adaptive Effector_memeory_CD8_T_cell immunity Gamma_delta_T_cell Immature__B_cell Regulatory_T_cell T_follicular_helper_cell Type_1_T_helper_cell Type_17_T_helper_cell Type_2_T_helper_cell Activated_dendritic_cell Immature_dendritic_cell CD56bright_natural_killer_cell Innate Macrophage immunity MDSC Monocyte Natural_killer_cell Natural_killer_T_cell Plasmacytoid_dendritic_cell BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC Normalized enrichment score (NES) –4 4 b Predictor variables (purity, age, and stage adjusted) Global methylation level Mutation burden Aneuploidy level Tumor type ALLOGRAFT_REJECTION COAGULATION COMPLEMENT Hallmark IL2_STAT5_SIGNALING immune IL6_JAK_STAT3_SIGNALING gene sets INFLAMMATORY_RESPONSE NTERFERON_ALPHA_RESPONSE INTERFERON_GAMMA_RESPONSE TNFA_SIGNALING_VIA_NFKB MAJOR_HISTOCOMPATIBILITY_COMPLEX ANTIGEN_PROCESSING_AND_PRESENTATION CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION E2F_TARGETS Hallmark G2M_CHECKPOINT proliferation MITOTIC_SPINDLE MYC_TARGETS_V1 gene sets MYC_TARGETS_V2 P53_PATHWAY BLCA BRCA CESC CRC ESCA HNSC KIRC KIRP LGG LIHC LUAD LUSC PAAD PCPG PRAD SARC SKCM STAD TGCT THCA UCEC NES –6 6 Fig. 2 Genomic methylation loss correlates with immune evasion signatures. a Heatmap of Gene Set Enrichment Analysis (GSEA)48 normalized enrichment scores (NESs) for gene sets representing various types of innate and adaptive immune cells (rows). For each gene per cancer type, a linear regression model was fit using mRNA expression level as response variable and global methylation, mutation burden, aneuploidy, tumour purity, age, and tumour stage as predictors. For each of three predictors (global methylation level, mutation burden, aneuploidy level), GSEA was performed on genes with significant regression coefficients. Cells with significant NES (FDR <0.25) are colour-scaled. b Heatmap of GSEA NESs for hallmark immune and proliferation gene sets and genes involved in antigen processing and presentation, MHC, or cytokine-cytokine receptor interaction expression levels of ERVs and L1s from the tumour samples. The long-range tumour demethylation was discovered in colon31, correlation of their expression levels with the indicators of breast32, and brain33 cancers. A recent study showed PMD cytotoxic immune activity was not positive but negative in general demethylation is a common feature of diverse cancer type16. Such (Supplementary Fig. 5), implying that IFN-α/β silencing overrides long-range demethylation in cancer is accompanied with gene the immune-stimulatory effects of ERV/L1 expression by genomic silencing programs. Genes within PMDs in differentiated cells are demethylation. under-expressed30. Similarly, genes in PMDs in various types of cancers are largely silenced by the formation of repressive chro- matin structures or via CGI hypermethylation31–35. Repression of immune genes in partial methylation domains. PMDs were characterized by and defined based on the high Methylation loss in late-replicating regions engages the formation variability of solo-WCGW methylation levels across samples16. of heterochromatic structure termed partial methylation domains Our inspection of the methylation variability and replication (PMDs) as opposed to highly methylated domains (HMDs)16. timing of various PMDs led to three distinct subclasses (Fig. 4a). PMDs were first discovered as contiguous regions with lower In accordance with a previous report36, the properties of PMDs levels of CpG methylation in differentiated cells30. PMD-like were associated with their genomic length with shorter PMDs 4 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 ARTICLE a Late-replicated genes (n = 885) BLCA BRCA CESC CRC ESCA HNSC LIHC LUAD LUSC PAAD PRAD STAD TGCT 1.0 mRNA expression level (average z-score) 0.5 0.0 −0.5 −1.0 P ** * * *** * ** *** ** ** ** *** *** ** ***P < 1.0 × 10–30 **P < 1.0 × 10–20 –10 Global L1 methylation level Low High *P < 1.0 × 10 b BLCA BRCA CRC HNSC KIRC KIRP LIHC LUAD LUSC PAAD PRAD SKCM STAD TGCT UCEC # of hyper-methylated CGIs 160 400 80 80 200 at late-replicated genes 300 150 150 150 300 200 100 100 20 80 200 40 100 200 40 100 50 50 100 50 50 50 10 100 P ** *** *** *** ** ** * *** *** *** *** *** *** *** *** Global L1 methylation level Low High ***P < 1.0 × 10–3 **P < 1.0 × 10–2 *P < 5.0 × 10–2 c d Cell cycle pathways Immune pathways (n = 108) (n = 130) Early-replicated genes (n = 1747) 20% 20 P = 3.7 × 10–9 25 Fraction of enriched pathways BLCA BRCA ESCA HNSC KIRP LUAD LUSC PAAD STAD UCEC 1.0 mRNA expression level (average z-score) 0.5 10% 0.0 6 −0.5 0 P ** *** ** *** * *** * * ** 0% ** Early replication Late replication –30 –20 –10 Global L1 methylation level Low High ***P < 1.0 × 10 **P < 1.0 × 10 *P < 1.0 × 10 in cancer in cancer e Enriched pathways in late-replicating regions (n = 6; FDR < 0.1) REGULATION_OF_IFNA_SIGNALING CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION NES 0 0 Enrichment score (ES) TRAF6_MEDIATED_IRF7_ACTIVATION CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION −3.0 −2.0 ES REGULATION_OF_IFNA_SIGNALING NES = –2.17 NES = –2.18 –3 JAK_STAT_SIGNALING_PATHWAY P = 9.0 × 10 P=0 −6.0 FDR = 0.01 −4.0 FDR = 0.02 RIG_I_MDA5_MEDIATED_INDUCTION_OF_IFN_ALPHA_BETA_PATHWAYS ADP_SIGNALLING_THROUGH_P2RY12 Earlier replication Later replication Earlier replication Later replication in cancer in cancer in cancer in cancer Fig. 3 Characterization of genes in late-replicating regions. a Comparison of the expression levels of late-replicated genes between tumours with low and high global methylation. Tumour types for which the paired t-test P < 1 × 10−10 are shown. b Comparison of the number of hypermethylated CpG island promoters of late-replicating genes between tumours with low and high global methylation. Tumour types with P < 0.05 (two-sided Mann–Whitney U test) are shown. In the boxplots, the centre line, bounds of box, and whiskers represent the 50th, 25th and 75th, and 5th, and 95th percentiles, respectively. c Comparison of the expression levels of early-replicated genes between tumours with low and high global methylation. Tumour types for which the paired t-test P < 1 × 10−10 are shown. d Fraction of enriched cell-cycle and immune pathways according to replication timing. The number of the pathways showing significant enrichment (FDR <0.25) by GSEA is indicated above the bars. e Significantly enriched pathways for genes in late-replicating regions and the two representative GSEA plots. Genes that are specifically expressed in the immune system were excluded characterized by earlier replication timing (Fig. 4a, b). Strikingly, Global methylation predicts responses to immunotherapy. To immunomodulatory pathway genes involved in antigen proces- test whether global methylation alterations affect the clinical sing and presentation, cytokine-cytokine receptor interaction, and benefit of immunotherapy, we generated methylome and exome JAK-STAT signaling pathway were concentrated in the short data for 60 samples in an anti-PD-1/PD-L1 cohort in lung cancer PMDs (Fig. 4c, d and Supplementary Table 2). The INF-α family collected from Samsung Medical Center (SMC) (Supplementary genes were in the short PMDs (Fig. 4e). Also, 8 HLA genes, Table 3). Also, we employed an additional anti-PD-1 lung cancer including HLA-DQA1, HLA-DRA, and HLA-DRB1, were located cohort composed of 81 methylomes and 22 exomes from Bellvitge within the short PMDs. Consistent with the late-replicating Biomedical Research Institute (IDIBELL)37. For validation, we regions (Fig. 3a, b), the short PMDs were accompanied with gene utilized data from 40 TCGA melanoma patients who received repression (Fig. 4f) and CGI hypermethylation (Fig. 4g) in immunotherapies. The summary of the three cohorts is provided demethylated tumours. Hypermethylated CGIs are most abun- in Fig. 5a. dant within 150 kb of PMD boundaries31. The enriched immune The samples from the combined lung cancer cohort were genes (Fig. 4c, d) were significantly concentrated near PMD divided into global low versus high methylation groups according boundaries with the average distance of 143 kb (Fig. 4h), to the L1 methylation levels. The global low methylation group suggesting that these genes are particularly prone to promoter exhibited decreased genomic (open sea/shelf) methylation and methylation. increased CGI/shore methylation (Fig. 5b, c). In agreement with NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 a b P = 1.7 × 10–64 d h Density 50 CELL_CYCLE 2e−06 6e−06 P = 1.5 × 10–34 P = 1.5 × 10–30 ANTIGEN_PROCESSING_AND_PRESENTATION P = 2.96 × 10–7 Average distance to the nearest HMD 150 kb Immune genes Long (L) PMD Short (S) PMD Intermediate (I) PMD 40 CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION (P = 0) 7.5 Mb P = 2.96 × 10–11 Domain size JAK_STAT_SIGNALING_PATHWAY Density 30 P = 2.09 × 10–4 250 kb Methylation 20 variability 2.5 Mb 10 350 kb Replication timing 0 0.05 0.10 0.15 0.20 Background S I L 450 kb distribution Replication timing Early Late Methylation variability Methylation variability Low High PMD class HMD S-PMD I-PMD L-PMD c e # of genes Class chr9p22 Enriched pathways in S-PMDs (%) S I L 21,050 kb 21,250 kb 21,450 kb AUTOIMMUNE_THYROID_DISEASE 22 (69%) ANTIGEN_PROCESSING_AND_PRESENTATION 28 (49%) Methylation HMD Short PMD domain CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION 47 (25%) GRAFT_VERSUS_HOST_DISEASE 14 (61%) IFNB1 IFNA4 IFNA10 IFNA17 IFNA5 IFNA2 IFNA1 REGULATION_OF_AUTOPHAGY 15 (56%) NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY 27 (28%) Genes INTESTINAL_IMMUNE_NETWORK_FOR_IGA_PRODUCTION 12 (41%) IFNA21 IFNA7 IFNA16 IFNA14 IFNA6 IFNA8 ASTHMA 9 (53%) RIG_I_LIKE_RECEPTOR_SIGNALING_PATHWAY 16 (30%) Early Normal cell lines (n n = 5) CYTOSOLIC_DNA_SENSING_PATHWAY 15 (31%) ALLOGRAFT_REJECTION 9 (45%) Replication timing JAK_STAT_SIGNALING_PATHWAY 24 (22%) TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY 18 (23%) Late Cancer cell lines (n n = 10) FDR < 0.01 f BLCA BRCA CRC KIRC LGG LIHC LUAD LUSC PAAD PCPG PRAD SKCM STAD TGCT THCA UCEC mRNA expression level 0.5 (average z-score) in short PMDs 0 –0.5 P * *** *** * * *** * * ** * ** *** * * ** * Global L1 methylation level Low High ***P < 1.0 × 10–30 **P < 1.0 × 10–20 *P < 1.0 × 10–10 g BLCA BRCA CESC CRC ESCA HNSC KIRP LGG LUAD LUSC PAAD PRAD SKCM STAD TGCT UCEC 200 125 150 200 150 200 150 100 100 80 150 150 80 200 120 # of hyper-methylated CGIs 120 in short PMDs 75 60 100 100 100 50 50 100 50 50 40 50 40 50 50 25 40 40 P *** *** *** *** *** *** * *** * *** *** * *** *** *** *** Global L1 methylation level Low High ***P < 1.0 × 10–3 **P < 1.0 × 10–2 *P < 5.0 × 10–2 Fig. 4 Characterization of genes in partially methylated domains. a Identification of PMD subclasses through hierarchical clustering on methylation variability and replication timing. b Comparison of PMD size between the identified PMD subclasses. c Enrichment of genes in immune-related pathways in the short PMDs. d Distribution of domain methylation variability for genes in cell cycle and immune pathways. The two-sample Kolmogorov-Smirnov test was used to assess deviation from the distribution of the cell cycle genes. e IFN-α genes in a short PMD with late replication timing. The mean and standard error of the weighted average signals of replication timing in normal cells and cancer cells are shown. f Comparison of the mRNA expression level of genes in the short PMDs between tumour samples with low and high global methylation. Tumour types for which the paired t-test P < 1.0 × 10−10 are shown. g Comparison of the number of hyper-methylated promoter CGIs in the short PMDs between tumour samples with low and high global methylation. Tumour types with P < 5.0 × 10−2 (two-sided Mann–Whitney U test) are shown. h Concentration of immune genes near PMD boundaries. The average distance of the immune-related genes (from b) to the nearest HMDs is marked by an arrow. The statistical significance of the observed average distance was assessed based on a null distribution generated by using random PMD genes our pan-cancer molecular data analyses, the global low methyla- Fig. 5g) resulted in HR = 0.57 (log rank test, P = 4.0 × 10−2) and tion samples showed high mutation burden and aneuploidy as HR = 0.52 (log rank test, P = 5.0 × 10−2), respectively. Whereas well as CGI hypermethylation in the short PMDs (Fig. 5d). the P values decreased as sample size grew, the effect size (HR) Transcriptome data of the SMC cohort and TCGA cohort showed remained a similar level. that genes involved in the MHC and cytokine-cytokine receptor We next compared the effect of global methylation and interaction were significantly enriched for repression in globally mutation load on the clinical response. For the combined lung demethylated tumours of both cohorts (Supplementary Fig. 6 and cancer cohort, methylome and matched exome data were Supplementary Data 3). available for 82 samples in total (Fig. 5a). In contrast to the High mutation load is associated with the clinical benefit of global L1 methylation level (log rank test, HR = 0.46 and P = checkpoint blockade2–5. However, when we stratified patient 7.0 × 10−3), mutation burden failed to show significant explana- samples according to the L1 methylation levels, the global low tory power (log rank test, HR = 0.77 and P = 0.3) (left and methylation group showed poor prognosis despite high mutation middle of Fig. 5h). Multiple regression with survival or clinical load. In the combined lung cancer cohort (n = 141), the hazard benefit as the response variable demonstrated significant effects ratio (HR) was 0.56 (log rank test, P = 7.0 × 10−3) (Fig. 5e). The by the global L1 methylation level but not mutation burden at IDIBELL cohort (n = 81) (Fig. 5f) and SMC cohort (left of P = 0.05 (right of Fig. 5h). We repeated the univariate 6 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 ARTICLE a b Methylation z-score level Global L1 methylation level –2 2 Low (n = 70) High (n=71) Number of samples Cancer Methylation Mutation Aneuploidy type Cohort level burden level SMC 60 60 60 Lung IDIBELL cancer 81 22 NA Combined 141 82 60 Melanoma TCGA 40 40 40 c Global L1 methylation d Low High P = 2.0 × 10–3 P = 4.3 × 10–9 P = 7.2 × 10–10 1e+09 # of hyper-methylated CGIs in short PMDs P=0 P=0 1 3000 200 Methylation level (β value) Aneuploidy level Mutation burden 1e+08 150 300 0.5 100 1e+07 30 50 0 CGI Open Sea Low High Low High Low High Probe 2 kb 2 kb 2 kb 2 kb (n = 27,497) (n = 199,030) location (n = 82) (n = 60) (n = 141) Open Sea Shelf Shore CGI Shore Shelf Open Sea Combined cohort (n = 141) IDIBELL cohort (n = 81) e Combined cohort (n = 141) f IDIBELL cohort (n = 81) g SMC cohort with methylation, mutation, and CNV data (n = 60) Global L1 methylation Global L1 methylation Global L1 methylation Mutation burden Aneuploidy level 1.0 1.0 1.0 1.0 1.0 HR = 0.52 HR = 0.64 HR = 1.56 HR = 0.56 P = 5.0 × 10–2 Progression free survival Progression free survival Progression free survival Progression free survival Progression free survival HR = 0.57 P = 0.2 P = 0.2 0.8 0.8 0.8 0.8 P = 7.0 × 10 –3 0.8 P = 4.0 × 10–2 0.6 0.6 0.6 0.6 0.6 0.4 0.4 0.4 0.4 0.4 High High Low 0.2 0.2 High 0.2 0.2 0.2 High Low High Low Low Low 0.0 0.0 0.0 0.0 0.0 0 10 20 30 40 0 10 20 30 40 0 5 10 15 20 25 0 5 10 15 20 25 0 5 10 15 20 25 Time (months) Time (months) Time (months) Time (months) Time (months) h Combined cohort with methylation and mutation data (n = 82) i IDIBELL cohort with methylation and mutation data (n = 22) Global L1 methylation Mutation burden 1.0 Global L1 methylation 1.0 Mutation burden Global L1 methylation 1.0 HR = 0.27 1.0 HR = 0.97 Mutaiton burden P = 2.0 × 10–2 Progression free survival Progression free survival Progression free survival HR = 0.46 HR = 0.77 Progression free survival P = 0.94 0.8 0.8 0.8 0.8 P = 7.0 × 10–3 P = 0.3 0.6 1.5 0.6 0.6 P = 0.05 0.6 -log10 (p value) 0.4 0.4 0.4 0.4 High High High 0.2 0.2 0.2 0.2 Low Low 0.5 Low Low High 0.0 0.0 0.0 0.0 0 5 10 15 20 25 0 5 10 15 20 PFS DCB/NDB 0 5 10 15 20 25 0 5 10 15 20 25 Time (months) Time (months) Time (months) Time (months) Fig. 5 Genomic demethylation adversely affects the clinical benefit of checkpoint blockade. a Summary of immunotherapy cohort samples with available data. b Unsupervised hierarchical clustering of the DNA methylation profiles of the SMC lung cancer cohort samples. The methylation probes (row) were clustered and samples (column) were divided into two groups according to the median of the global methylation level. The heatmap shows beta values for the most differentially methylated loci (the highest 1%) between tumour samples with low and high global methylation. The methylation probes were categorized by the relative distance from CGIs (CGI, shore, shelf, and open sea). c Differential CGI/open-sea methylation between the global low and high groups. The differentially methylated loci between the two groups (FDR <0.05) were selected, and per-locus beta values were averaged for each group. The paired t-test was used to test the statistical significance. d Comparison of mutation burden, the aneuploidy level, and the number of hyper-methylated promoter CGIs in the short PMDs between tumour samples with low and high global methylation. e, f Survival analysis using tumour samples with methylation data from the e combined or f IDIBELL cohort. Patients were stratified by the global L1 methylation level. The log rank test was used to compare survival curves estimated by the Kaplan–Meier method. g Survival analysis using the SMC tumour samples with methylation level, mutation burden, and aneuploidy level. Patients were stratified by the global L1 methylation level (left), mutation burden (middle), and the aneuploidy level (right). h Survival analysis using tumour samples with methylation and mutation data from the combined cohort. Patients were stratified by the global L1 methylation level (left) or mutation burden (middle). To estimate the relative contribution of mutation burden and global methylation to patient survival and clinical benefit, the multivariate Cox proportional hazards model (for PFS: progression free survival) or multivariate logistic regression model (for DCB: durable clinical benefit and NDB: no durable benefit). i Survival analysis using the IDIBELL tumour samples with methylation and mutation data. Patients were stratified by the global L1 methylation level (left) or mutation burden (right) NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 comparisons for the two lung cancer cohorts separately. Both the for their own elimination by the immune system as a means for SMC cohort (n = 60) and IDIBELL cohort (n = 22) resulted in cancer cell immunosurveillance. significant stratification by global methylation but not by mutation burden (left and middle of Fig. 5g and i). Discussion The global L1 methylation level negatively correlated with In this work, we propose that DNA methylation aberration is an aneuploidy (middle of Fig. 5d). A pan-cancer analysis discovered important determinant of the tumour response to host immune the association of aneuploidy with signatures of immune activity, and can provide a mechanism by which rapidly dividing evasion10. Taken together, both low methylation and high and highly mutated tumours escape immune reaction and resist aneuploidy are expected to decrease tumour immunity and immunotherapy. The key mechanism seems to be the formation undermine the clinical benefit of immunotherapy. However, of heterochromatin, which is coupled with progressive domain- when we examined the SMC lung cancer cohort for which level methylation loss. An open question is what dictates these aneuploidy data were available (Fig. 5a), only global methylation epigenetic changes in particular regions, such as the MHC locus. but not aneuploidy showed significant correlations with poor One possibility is that these changes at the particular loci are clinical responses (left and right of Fig. 5g). selectively favoured during cancer evolution because they provide All the above analyses were repeated for the TCGA melanoma immune evasion mechanisms and increase fitness of tumour cells. cohort (n = 40) (Fig. 5a). The samples were first divided into It is also possible that domain demethylation of immune-related global low versus high methylation groups according to the L1 regions is a more inherent chromatin tendency than a con- methylation levels. The global low methylation group in sequence of selection. In any case, our results suggest that mitotic melanoma recapitulated CGI promoter hypermethylation in the cell division causes genetic and epigenetic alterations that exert short PMDs (Supplementary Fig. 7A) and poor prognosis in opposing effects on tumour immunity by increasing neoantigens response to immunotherapies (log rank test, HR = 0.48 and P = and inhibiting immune gene expression, respectively. Cell divi- 3.0 × 10−2) (Supplementary Fig. 7B). Mutation burden and sion also increases focal and chromosomal copy number changes. aneuploidy level both failed to explain the clinical benefit Our data show that the particular association of chromosomal (Supplementary Fig. 7C, D). copy number changes with low antitumour immune activity can be explained by global methylation loss. There are multiple studies that reported antitumour immunity Global demethylation rules out the effect of aneuploidy. Our augmented by CDK4/6 inhibition and synergistic effects of the results on aneuploidy from the clinical data contradict the previous cell cycle inhibitors and checkpoint blockade44–46. Our results findings on the association of aneuploidy and immune evasion suggest that cell cycle inhibition may bring about opposing effects signatures10. Thus, we examined the possibility that global by suppressing genetic alterations that facilitate neoantigen for- methylation is able to explain this association. Global demethyla- mation and at the same time, preventing immune evasion pro- tion in cancer promotes chromosomal instability17–20. DNA moted by epigenetic alterations that repress immunomodulatory hypomethylation-related instability is mainly of chromosomal pathway genes. Hence, the reported effects of cell cycle inhibition nature and involves large-scale alterations leading to aneuploidy suggest that the benefits achieved by epigenetic influences may be rather than widespread amplifications or deletions21–23. As in ICF greater than the adverse effects caused by suppressing neoantigen syndrome (for immunodeficiency, centromere instability, and formation. facial anomalies), failure of methylation maintenance in pericen- The repression signatures for IFN-α/β signaling draw parti- tromeric sequences can cause erroneous chromosomal segregation cular attention, given that this pathway is supposed to stimulate in cancer38–42. immune responses against dsRNAs induced by genomic deme- Indeed, global demethylation significantly correlated with thylation. Recent studies have shown that DNA methylation SCNAs across different tumour types (Fig. 6a). We determined inhibitors induce dsRNA expression and stimulate antitumour the magnitude of chromosomal SCNAs (cSCNAs) by combining immune activity through the IFN-α/β response activated by the the chromosome SCNA and arm SCNA levels that were viral defence pathway27–29. Based on these results, combining previously calculated10, and compared this result with that of epigenetic therapy and immunotherapy has been suggested47,48. focal SCNAs (fSCNAs). The correlation was stronger with According to our results, tumours with global methylation loss cSCNAs than with fSCNAs (Fig. 6b and Supplementary Fig. 8). tend to resist immunotherapy alone and may particularly benefit We used partial correlations to estimate the extent to which the from this combined treatment approach. However, intrinsically global L1 methylation level correlates with cSCNAs (or fSCNAs) de-suppressed dsRNAs may fail to boost antitumour immunity when controlling for fSCNAs (or cSCNAs). In most cases, global because of inactivated IFN-α/β signaling. We indeed observed hypomethylation was associated with cSCNAs independently of that ERV/LINE expression does not increase the immune sig- fSCNAs (Fig. 6c). In contrast, the correlation with fSCNAs natures, which, in contrast, are reduced probably reflecting IFN- disappeared when cSCNAs were controlled for (Fig. 6c). α/β inactivation. Therefore, different action mechanisms of Importantly, the multiple regression analysis of the immune epigenetic therapy are required when targeting these tumours. signature scores revealed markedly higher explanatory power for Specifically, it needs to be tested whether methylation inhibitors the global L1 methylation levels than cSCNA levels (Fig. 6d). We or other epigenetic modulators are capable of restoring the IFN- also performed the partial correlation analyses for the immune α/β response and other immunomodulatory pathways by signature score, global L1 methylation level, and cSCNA level. diminishing CGI methylation or loosening heterochromatin Overall, the positive correlation between the immune signature structure in these intrinsically demethylated tumours. Our study score and global L1 methylation level maintained when the sheds light on the combination of epigenetic modulation and cSCNA level was controlled for (Fig. 6e), except for one tumour checkpoint blockade as a potential precision immunotherapy type that was previously reported as an exception regarding the regimen. role of aneuploidy10. These results indicate that the immune avoidance signatures of highly aneuploid tumours are associated with genomic demethylation. Indeed, a recent molecular Methods TCGA molecular and clinical data. The batch-corrected and normalized DNA mechanism study43 contradicted the previous report10 by methylation data based on Infinium Methylation 450k technology), together with suggesting that aneuploid cells generate pro-inflammatory signals mRNA expression and gene mutation data, generated by the PanCancer Atlas 8 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 ARTICLE a BRCA CRC GBM HNSC KIRC LGG LUAD LUSC SKCM STAD UCEC P = 3.7 × 10–12 P = 1.2 × 10–10 P = 3.1 × 10–2 P = 1.6 × 10–19 P = 3.2 × 10–4 P = 3.9 × 10–7 P = 3.1 × 10–21 P = 2.8 × 10–16 P = 5.7 × 10–11 P = 5.5 × 10–19 P = 2.5 × 10–8 2 SCNA level 0 −2 Low High Low High Low High Low High Low High Low High Low High Low High Low High Low High Low High n = 334 n = 194 n = 68 n = 194 n = 164 n = 224 n = 198 n = 124 n = 136 n = 144 n = 206 Global DNA methylation level (lowest 30% vs highest 30%) b d A EC M SC AD SC C AD BM C G C C R R BRCA LG BR C N LU LU ST SK KI G C U H CRC Chromosomal SCNA level 2 2 GBM ** ** ** ** ** ** ** ** ** Focal SCNA level Methylation 7.27 3.48 1.78 5.13 1.71 –3.45 3.53 4.36 4.63 4. 11 4.25 HNSC KIRC * * * * 0 0 ** LGG cSCNA –3.08 –3.33 0.97 –1.04 –0.42 2.99 –3.09 –1.04 –2.91 –2.49 –2.21 LUAD * −2 −2 LUSC fSCNA 3.24 -2.14 0.89 –1.70 0.86 –0.59 1.33 0.10 2.14 0.05 –0.16 r = –0.43 r = –0.33 SKCM –153 P = 1.3 × 10 n = 3342 P = 1.7 × 10–87 n = 3342 STAD Z **P < 1.0 × 10–3 *P < 1.0 × 10–2 −2 0 2 −2 0 2 UCEC 7.5 0 –7.5 Global methylation level Global methylation level c e BRCA CRC GBM HNSC KIRC LGG LUAD LUSC SKCM STAD UCEC BRCA CRC GBM HNSC KIRC LGG LUAD LUSC SKCM STAD UCEC 0 0 0 0 0 0 0 0 0 0 0.2 0.5 0 0.5 0.5 0.2 0.2 0.2 0.2 0.2 0.2 0.2 r 0 0 0 0 0 0 0 0 0 0 0 −0.2 −0.2 –0.3 –0.15 –0.15 –0.15 –0.3 –0.4 –0.4 –0.3 –0.1 –0.2 –0.2 –0.2 –0.2 –0.2 –0.2 –0.5 –0.5 –0.2 –0.2 –0.5 P *** *** * *** ** *** *** *** *** *** *** *** ** *** ** * *** *** *** * *** *** * *** * *** *** *** *** *** *** *** *** * cSCNA level fSCNA level ***P < 1.0 × 10–3 **P < 1.0 × 10–2 *P < 5.0 × 10–2 Methylation level cSCNA level ***P < 1.0 × 10–3 **P < 1.0 × 10–2 *P < 5.0 × 10–2 Fig. 6 Aneuploidy indicates global DNA hypomethylation. a Comparison of SCNA levels between tumours with low and high global methylation. Shown here is the combination of cSCNAs and fSCNAs. P values from the two-sided Mann–Whitney U test are shown. b Association of global DNA methylation with cSCNA and fSCNA levels. The values were standardized per tumour type. Spearman’s rank correlation coefficient and its P value are shown. c Partial correlation analyses comparing cSCNAs (red) and fSCNAs (orange) for their correlation with global methylation independently of one another. The Spearman correlation was used. d For each tumour type, samples with a low (<30th percentile) and high (>70th percentile) immune signature score were compared by multivariate logistic regression with the standardized global methylation, cSCNA, and fSCNA levels, and tumour purity as predictors. e Partial correlation analyses comparing global methylation (blue) and cSCNAs (red) for their correlation with the immune signature score independently of one another. The Spearman correlation was used. (a–e) These analyses were confined to 11 tumour types for which the cSCNA and fSCNA data were available consortium were obtained from the publication page (https://gdc.cancer.gov/about- Supplementary Table 5) from the GDC legacy archive (https://portal.gdc.cancer. data/publications/pancanatlas). TCGA aliquot barcodes flagged as DO NOT USE gov/legacy-archive). We first selected evolutionarily young LINE-1 repeat elements in the Merged Sample Quality Annotation file were discarded. We selected cancer containing at least three different CpG sites that were covered by at least 10 aligned types for which there were >100 patient samples with all the molecular data and reads and then computed the averaged beta value for each repeat element. The age information. The selected tumour types encompassed 6968 samples across 21 number of the selected repeat elements for each sample ranged from 67,732 to types (Supplementary Table 4), including bladder urothelial carcinoma (BLCA), 71,958 (Supplementary Table 5). We obtained the mean of the average beta values breast adenocarcinoma (BRCA), cervical squamous cell carcinoma and endo- of the repeat elements to represent the global methylation level of each sample. We cervical adenocarcinoma (CESC), colorectal carcinoma (CRC), esophageal carci- performed an additional validation by using the LINE-1 pyrosequencing dataset of noma (ESCA), head and neck squamous cell carcinoma (HNSC), kidney renal clear 15 samples for which methylation array data of the same platform were available52. cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), lower grade We compared our measures based on the selected array probes with the bisulfite glioma (LGG), liver hepatocellular carcinoma (LIHC), lung adenocarcinoma sequencing and pyrosequencing measures (Supplementary Fig. 2). (LUAD), lung squamous cell carcinoma (LUSC), pancreatic adenocarcinoma (PAAD), pheochromocytoma and paraganglioma (PCPG), prostate adenocarci- noma (PRAD), sarcoma (SARC), skin cutaneous melanoma (SKCM), stomach Linear regression modelling. For each gene per cancer type, a linear regression adenocarcinoma (STAD), testicular germ cell tumours (TGCT), thyroid carcinoma model was fit using mRNA expression level as the response variable, and global (THCA), and uterine corpus endometrial carcinoma (UCEC). PMD methylation methylation, mutation burden, aneuploidy, tumour purity, age, and tumour stage levels were derived from Zhou et al.16 (https://zwdzwd.github.io/pmd). Aneuploidy as predictors. We included tumour stage in the model for tumour types for which level and tumour purity values were obtained from a Table of Taylor et al.49. The at least 100 patient samples with tumour stage information were available (n = 15; cSCNA and fSCNA levels from 11 cancer types were obtained from Davoli et al.10. Supplementary Table 4). The regression model with the following formula was For glioblastoma multiforme (GBM), microarray-based gene expression data were built using the lm function in R. obtained from the UCSC Xena public data hubs. mRNA expression of gene Y ~ β1 × global methylation level + β2 × mutation burden + β3 × aneuploidy level + β4 × tumour purity + β5 × age + β6 × tumour stage For each of three predictors (global methylation level, mutation burden, aneuploidy level) per tumour type, GSEA48 was performed on genes with Estimating global methylation levels. To determine global methylation levels, we significant regression coefficients (Benjamini and Hochberg FDR <0.05). Genes chose methylation probes for which at least 90% of sequences (≥45 bp) mapped to ranked by the z score from the linear regression model were used for input into the the young subfamilies of LINE-1 (L1HS and L1PA). Probe mapping information preranked module of the GSEA software48 with the immune and proliferation gene and LINE-1 family annotation were obtained from GPL1630450 and the Repeat- sets (see below). GSEA was run with default settings, except for the minimum Masker of the UCSC genome browser51, respectively. We averaged the beta values number of gene sets, which was set to 10. of the chosen probes in each tumour sample, and used the average beta as an estimate of the global methylation level (Supplementary Data 1). For validation, we obtained the whole-genome bisulfite sequencing data (level 3 beta values) of Identification of genes with differential replication timing. Repli-Seq mea- samples for which the methylation array data were available (n = 18; surements (wavelet-smoothed signal) of ENCODE 5 normal (HUVEC, IMR90, NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 NHEK, BJ_1 and BJ_2) and 10 cancer cell lines (MCF-7, SK-N-SH, HepG2, HeLa, responder. Progression-free survival (PFS) was calculated from the start date of A549, G401, LNCaP, T47D, H460, and Caki-2) were downloaded from the UCSC therapy to the date of progression or death, whichever is earlier. Patients were Genome Browser53 and ENCODE project portal54. Average wavelet-smoothed censored at the date of the last follow-up for PFS if they were not progressed and signal values in each 5-kb window were scaled and quantile-normalized. Windows alive. We complied with all relevant ethical regulations for work with human par- harbouring missing values in any of the cell lines were excluded. We performed the ticipants. Informed consent was obtained. This study was approved from the insti- student’s t-test for each 5-kb window to assess the replication timing difference tutional review board at Samsung Medical Center (2018-03-130 and 2013-10-112). between the normal and cancer cells. We then assigned the P value and t statistic to each gene. For genes spanning multiple replication timing windows, we assigned a combined P value using the Fisher’s method and average t statistic. We excluded Whole-exome, transcriptome, and methylome data for the SMC cohort. genes located on sex chromosomes or specifically expressed in the immune system Tumour samples were obtained before anti-PD1/PD-L1 treatment, and then were (see below). Genes with Bonferroni-adjusted P value < 0.05 were defined as early- embedded in paraffin after formalin fixation or kept fresh. DNA was prepared or late-replicated genes in cancer. Genes ranked by the t statistic were used for using AllPrep DNA/RNA Mini Kit (Qiagen, 80204), AllPrep DNA/RNA Micro Kit input into the preranked module of the GSEA software48 with the canonical (Qiagen, 80284), or QIAamp DNA FFPE Tissue Kit (Qiagen, 56404) for library REACTOME and KEGG pathways from MSigDB55. preparation for whole exome sequencing. Library preparation was performed by using SureSelectXT Human All Exon V5 (Agilent, 5190–6209) according to the instructions61. Briefly, 200–300 ng of tumour and normal genomic DNA was Analysis of genes with differential replication timing. To investigate the sheared, and 150–200 bp of the sheared DNA fragments were further processed for expression differences of early- and late-replicated genes between tumour samples end-repairing, phosphorylation, and ligation to adaptors. Ligated DNA was with high and low global methylation levels, we partitioned tumour samples into hybridized using whole-exome baits from SureSelectXT Human All Exon V5. The the low (<30th percentile) and high (>70th percentile) global methylation groups for libraries were quantified by Qubit and 2200 Tapestation, and sequenced on an each cancer type. After a z-score normalization of mRNA expression data per gene Illumina HiSeq 2500 platform with 2 × 100 bp paired ends. Target coverage for per cancer type, we computed the average expression level of the genes for each normal samples was × 50 and tumour sample was × 100. group and compared the groups. To calculate the fraction of enriched cell-cycle The sequencing reads were aligned to the human reference genome (hg19) with and immune pathways (Fig. 3c), genes ranked by the t statistic (replication timing BWA mem module (v0.7.12)60 with default parameters. PCR duplicate reads were difference between normal and cancer cells) were used for input into the preranked marked using Picard62. We used Strelka263 to call somatic variants and selected module of the GSEA software with REACTOME pathways (provided at https:// single nucleotide variants (SNVs) and indels covered by at least ten and five reads reactome.org/) belonging to the Cell Cycle or Immune System category. Pathways in tumour, respectively. We further filtered out common germline variants present harbouring at least 10 genes were used for this analysis. in dbSNP 15064 and annotated somatic variants using ANNOVAR65. The list of filtered SNVs and indels is provided in Supplementary Data 5. Copy number variations (CNVs) were called using CNVkit66 with the circular binary Clustering of PMD and analysis of focal hypermethylation in PMD. By segmentation algorithm (Supplementary Data 6). Aneuploidy levels were derived employing previously defined locations of PMDs and HMDs per 100 kb16, we from the called CNVs. Specifically, we applied the defined threshold of ± 0.2 merged consecutive domains of the same type and retained those >300 kb in (average value of LUAD and LUSC) on the segment log2 ratio (tumour versus length. After assigning average methylation variability and replication timing from normal) to detect amplifications/deletions affecting at least 10% of a chromosome normal samples into each merged PMD, we performed hierarchical clustering on arm or 5% of a chromosome. The aneuploidy level was the sum of the absolute them. The enrichment of genes in particular pathways for clusters was computed segment log2 ratio, each weighted by its length10. using the binomial test. To estimate the number of hypermethylated CpG island RNA was extracted from same tumour tissue using Allprep DNA/RNA Mini promoter probes in short PMDs for each sample, we first calculated the mean and Kit (Qiagen, 80204). RNA was extracted from formalin fixed paraffin embedded standard deviations of CpG island probes in the merged HMDs and then counted (FFPE) using Rneasy FFPE kit (Qiagen, 73504). RNA was assessed for quality and CpG island promoter probes (annotated as TSS200 or TSS1500) in short PMDs for quantity using nanodrop 8000 UV-Vis spectrometer (NanoDrop Technologies Inc) which methylation level is greater than two standard deviations from the mean. and 4200 TapeStation Instrument (Aglient Technologies). RNA integrity number (RIN) of > = 5 were selected for further library preparation. In total 500 ng of RNA Analysis of proximity of immune genes to PMD boundaries. For a total of 77 from fresh tissues and 100 ng of RNA from FFPE were used for library preparation immune genes in the pathways enriched for short PMDs (n = 13; Fig. 4c), we using Truseq RNA library prep kit v2 (Illumina, RS-122-2001, Rs-122-2002) or calculated the average distance to their nearest HMD (original HMD defined by Truseq RNA access library prep kit (Illumina, RS-301-2001, RS-301-2002), Zhou et al.16) and then estimated its P value by generating a background dis- respectively. The library was generated according to the manufacturer’s tribution. We randomly picked 77 PMD genes and calculated the average distance instructions. RNA libraries were multiplexed and sequenced with 100 bp pair end of them to the nearest HMD. This procedure was repeated 10,000 times. reads on HiSeq2500 platform (Illumina). The RNA-seq reads were aligned to the human reference genome (hg19) with STAR67 and gene expression values were quantified using RSEM68. Genes ranked Collection and identification of marker gene sets. We obtained markers for by t-values obtained from comparing mRNA expression levels between tumours CD8 + T cells and proliferation (Fig. 1a, d) from Thorsson et al.56 and used single- with low (n = 14) and high global methylation level (n = 13) were used for input sample GSEA to estimate the activity of the markers. Gene sets for marking dif- into the preranked module of the GSEA software with KEGG pathways and the ferent types of immune cells and MHC (class I, class II, and non-classical) were MHC gene set. derived from Charoentong et al.57. Hallmark immune and proliferation gene sets Methylation assay was performed by following the instructions of Infinium were obtained from MSigDB55. Antigen presentation and cytokine signaling MethylationEPIC BeadChIP Kit (Illumina, WG-317-1002). Briefly, 500 ng genomic pathways were derived from canonical KEGG pathway in MSigDB55. The gene set DNA (gDNA) was used for bisulfite conversion using the EZ DNA methylation kit for the immune signature score was obtained from the aneuploidy study10. To filter (Zymo Research, D5001). The bisulfited gDNA was denatured and neutralized for out genes that are specifically expressed in the immune system, we used gene amplification, and was further processed for fragmentation. After fragmentation, expression data from the Illumina’s Human BodyMap 2.0 project (ftp://ftp.ncbi. DNA was eluted and resuspended in a hybridization buffer, and then hybridized nih.gov/gene/DATA/expression/Mammalia/Homo_sapiens/). Genes for which the onto the BeadChip. The BeadChip was prepared for staining and extension after average expression level in leukocytes and lymph nodes was five-fold higher than washing out unhybridized DNA, and it was imaged using the Illumina iScan System. that in the remaining tissues (n = 14) were considered as genes specifically The raw intensity files were then preprocessed into beta values using the expressed in the immune system (n = 1216; Supplementary Data 4). preprocessIllumina function in minfi69. The methylation data were treated as described in the Estimating global methylation levels section. The PMD levels of our Quantification of LINE-1 and ERV expression. To quantify the LINE-1 and ERV cohort samples were calculated based on the average of EPIC probes for Solo- expression levels, we aligned RNA-seq reads against the LINE-158 and ERV59 WCGW CpGs in common PMDs16 (provided at https://zwdzwd.github.io/pmd). sequence library, respectively, by using BWA60. We then normalized the mapped Redundant probes such as multi-hit probes by using the filter function of the read counts by the total number of aligned RNA-seq reads. Reads that mapped to ChAMP package70. We used MethylCIBERSORT71 and ESTIMATE72 to estimate both of the libraries or other repeat libraries58 (Alu and SVA) were excluded. The tumour purity and leukocyte fraction (Supplementary Fig. 4). We processed the raw normalized expression levels were standardized per tumour type for comparison. methylation intensity files of 81 lung cancer samples of the IDIBELL cohort37 with the same pipeline and merged them with the SMC cohort data using ComBat73. Clinical data of the SMC cohort. A total of 60 advanced non-small cell lung carcinoma patients who were treated with anti-PD-1/PD-L1 from 2014 to 2017 at Melanoma cohort data. Progression-free survival data for melanoma patients who Samsung Medical Center were enrolled for this study (Supplementary Table 3). The received immune checkpoint inhibitors (drug name labelled as Ipilimumab, Yer- clinical response was evaluated by the Response Evaluation Criteria in Solid voy, or Pembrolizumab; n = 15) were obtained from Ock et al.9. We included Tumours (RECIST) version 1.1 with a minimum 6-month follow-up. The response additional 25 patients that received other types of immunotherapy using drug data to immunotherapy was classified into durable clinical benefit (DCB, responder) or from the GDC legacy archive. We selected samples for which the therapy type non-durable benefit (NDB, non-responder)2. Partial response (PR) or stable disease (CDE_ID:2793530) column indicated immunotherapy while excluding samples (SD) that lasted more than 6 months was considered as DCB/responder. Progressive from patients that received multiple drugs. The molecular data for these samples disease (PD) or SD that lasted less than 6 months was considered as NDB/non- were obtained as described in the TCGA molecular and clinical data section. 10 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 ARTICLE Multivariate survival analysis. Global methylation and mutation burden were 23. Suzuki, K. et al. Global DNA demethylation in gastrointestinal cancer is age combined in a multivariable Cox proportional hazards model using the coxph dependent and precedes genomic damage. Cancer Cell 9, 199–207 (2006). function in R. The multivariable logistic regression model was used to assess the 24. Weisenberger, D. J. et al. Analysis of repetitive element DNA methylation by impact of global methylation and mutation burden on the objective response using MethyLight. Nucleic Acids Res. 33, 6823–6836 (2005). the glm function in R. 25. Yang, A. S. A simple method for estimating global DNA methylation using bisulfite PCR of repetitive DNA elements. Nucleic Acids Res. 32, 38e–38e Reporting summary. Further information on research design is available in (2004). the Nature Research Reporting Summary linked to this article. 26. Reis, E. S., Mastellos, D. C., Ricklin, D., Mantovani, A. & Lambris, J. D. Complement in cancer: untangling an intricate relationship. Nat. Rev. Immunol. 18, 5–18 (2018). Data availability 27. Chiappinelli, K. B. et al. Inhibiting DNA methylation causes an interferon The methylation chip and RNA-seq data for the samples of our lung cancer cohort are response in cancer via dsRNA including endogenous retroviruses. Cell 162, available at Gene Expression Omnibus under GSE119144 and GSE135222, respectively. 974–986 (2015). The raw data for the exome sequencing of our SMC cohort samples have been submitted 28. Roulois, D. et al. DNA-demethylating agents target colorectal cancer cells by to European Genome-phenome Archive under accession number EGAS00001003731. inducing viral mimicry by endogenous transcripts. Cell 162, 961–973 (2015). 29. Topper, M. J. et al. Epigenetic therapy ties MYC depletion to reversing Code availability immune evasion and treating lung cancer. Cell 171, 1284–1300 (2017). e21. Computer codes used in this study are provided as Supplementary Software 1. 30. Lister, R. et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature 462, 315–322 (2009). 31. Berman, B. P. et al. Regions of focal DNA hypermethylation and long-range Received: 24 October 2018 Accepted: 10 August 2019 hypomethylation in colorectal cancer coincide with nuclear lamina–associated domains. Nat. Genet. 44, 40–46 (2012). 32. Hon, G. C. et al. Global DNA hypomethylation coupled to repressive chromatin domain formation and gene silencing in breast cancer. Genome Res. 22, 246–258 (2012). 33. Hovestadt, V. et al. Decoding the regulatory landscape of medulloblastoma References using DNA methylation sequencing. Nature 510, 537–541 (2014). 1. Sharma, P. & Allison, J. P. The future of immune checkpoint therapy. Science 34. Timp, W. et al. Large hypomethylated blocks as a universal defining epigenetic 348, 56–61 (2014). alteration in human solid tumors. Genome Med. 6, 61 (2014). 2. Rizvi, N. A. et al. Mutational landscape determines sensitivity to PD-1 35. Brinkman, A. B. et al. Partially methylated domains are hypervariable in breast blockade in non-small cell lung cancer. Science 348, 124–128 (2015). cancer and fuel widespread CpG island hypermethylation. Nat. Commun. 10, 3. Schumacher, T. N. & Schreiber, R. D. Neoantigens in cancer immunotherapy. 1749 (2019). Science 348, 69–74 (2015). 36. Salhab, A. et al. A comprehensive analysis of 195 DNA methylomes reveals 4. Gubin, M. M. & Schreiber, R. D. The odds of immunotherapy success. Science shared and cell-specific features of partially methylated domains. Genome Biol. 350, 158–159 (2015). 19, 150 (2018). 5. Van Allen, E. M. et al. Genomic correlates of response to CTLA-4 blockade in 37. Davalos, V. et al. Epigenetic prediction of response to anti-PD-1 treatment in metastatic melanoma. Science 350, 207–211 (2015). non-small-cell lung cancer: a multicentre, retrospective analysis. Lancet Respir. 6. Tumeh, P. C. et al. PD-1 blockade induces responses by inhibiting adaptive Med. 6, 771–781 (2018). immune resistance. Nature 515, 568–571 (2014). 38. Costa, G., Barra, V., Lentini, L., Cilluffo, D. & Leonardo, A. Di. DNA 7. Herbst, R. S. et al. Predictive correlates of response to the anti-PD-L1 antibody demethylation caused by 5-Aza-2’-deoxycytidine induces mitotic alterations MPDL3280A in cancer patients. Nature 515, 563–567 (2014). and aneuploidy. Oncotarget 7, 3726–3739 (2016). 8. Roh, W. et al. Integrated molecular analysis of tumor biopsies on sequential 39. Tsuda, H., Takarabe, T., Kanai, Y., Fukutomi, T. & Hirohashi, S. Correlation CTLA-4 and PD-1 blockade reveals markers of response and resistance. Sci. of DNA hypomethylation at pericentromeric heterochromatin regions of Transl. Med. 9, eaah3560 (2017). chromosomes 16 and 1 with histological features and chromosomal 9. Ock, C. Y. et al. Genomic landscape associated with potential response to anti- abnormalities of human breast carcinomas. Am. J. Pathol. 161, 859–866 CTLA-4 treatment in cancers. Nat. Commun. 8, 1050 (2017). (2002). 10. Davoli, T., Uno, H., Wooten, E. C. & Elledge, S. J. Tumor aneuploidy 40. Wong, N. et al. Hypomethylation of chromosome 1 heterochromatin DNA correlates with markers of immune evasion and with reduced response to correlates with q-arm copy gain in human hepatocellular carcinoma. Am. J. immunotherapy. Science 355, 8399 (2017). Pathol. 159, 465–471 (2001). 11. McGranahan, N. et al. Clonal neoantigens elicit T cell immunoreactivity and 41. Qu, G. Z., Grundy, P. E., Narayan, A. & Ehrlich, M. Frequent sensitivity to immune checkpoint blockade. Science 351, 1463–1469 (2016). hypomethylation in Wilms tumors of pericentromeric DNA in chromosomes 12. Riaz, N. et al. Recurrent SERPINB3 and SERPINB4 mutations in patients who 1 and 16. Cancer Genet. Cytogenet. 109, 34–39 (1999). respond to anti-CTLA4 immunotherapy. Nat. Genet. 48, 1327–1330 (2016). 42. Xu, G.-L. et al. Chromosome instability and immunodeficiency syndrome 13. Gao, J. et al. Loss of IFN-γ pathway genes in tumor cells as a mechanism of caused by mutations in a DNA methyltransferase gene. Nature 402, 187–191 resistance to Anti-CTLA-4 Therapy. Cell 167, 397–404 (2016). (1999). 14. Shipony, Z. et al. Dynamic and static maintenance of epigenetic memory in 43. Santaguida, S. et al. Chromosome Mis-segregation generates cell-cycle- pluripotent and somatic cells. Nature 513, 115–119 (2014). arrested cells with complex karyotypes that are eliminated by the immune 15. Ehrlich, M. DNA methylation in cancer: too much, but also too little. system. Dev. Cell 41, 638–651 (2017). e5. Oncogene 21, 5400–5413 (2002). 44. Goel, S. et al. CDK4/6 inhibition triggers anti-tumour immunity. Nature 548, 16. Zhou, W. et al. DNA methylation loss in late-replicating domains is linked to 471–475 (2017). mitotic cell division. Nat. Genet. 50, 591–602 (2018). 45. Deng, J. et al. CDK4/6 inhibition augments antitumor immunity by enhancing 17. Jaenisch, R., Chen, R. Z., Pettersson, U., Beard, C. & Jackson-Grusby, L. DNA T-cell activation. Cancer Disco. 8, 216–233 (2018). hypomethylation leads to elevated mutation rates. Nature 395, 89–93 (1998). 46. Schaer, D. A. et al. The CDK4/6 inhibitor abemaciclib induces a T cell 18. Eden, A., Gaudet, F., Waghmare, A. & Jaenisch, R. Chromosomal instability inflamed tumor microenvironment and enhances the efficacy of PD-L1 and tumors promoted by DNA hypomethylation. Science 300, 455–455 checkpoint blockade. Cell Rep. 22, 2978–2994 (2018). (2003). 47. Chiappinelli, K. B., Zahnow, C. A., Ahuja, N. & Bylin, S. B. Combining 19. Gaudet, F. et al. Induction of tumors in mice by genomic hypomethylation. epigenetic and immunotherapy to combat cancer. Cancer Res. 76, 1683–1689 Science 300, 489–492 (2003). (2016). 20. Lee, S. T. & Wiemels, J. L. Genome-wide CpG island methylation and 48. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based intergenic demethylation propensities vary among different tumor sites. approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Nucleic Acids Res. 44, 1105–1117 (2016). Sci.USA 102, 15545–15550 (2005). 21. Karpf, A. R. & Matsui, S. I. Genetic disruption of cytosine DNA 49. Taylor, A. M. et al. Genomic and functional approaches to understanding methyltransferase enzymes induces chromosomal instability in human cancer cancer aneuploidy. Cancer Cell 33, 676–689 (2018). e3. cells. Cancer Res. 65, 8635–8639 (2005). 50. Price, M. E. et al. Additional annotation enhances potential for biologically- 22. Rodriguez, J. et al. Chromosomal instability correlates with genome-wide relevant analysis of the Illumina Infinium HumanMethylation450 BeadChip DNA demethylation in human primary colorectal cancers. Cancer Res. 66, array. Epigenetics Chromatin 6, 1–15 (2013). 8462–9468 (2006). 51. Rosenbloom, K. R. et al. The UCSC Genome Browser database: 2015 update. Nucleic Acids Res. 43, D670–D681 (2015). NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-12159-9 52. Bock, C. et al. Quantitative comparison of DNA methylation assays for 73. Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray biomarker development and clinical applications. Nat. Biotechnol. 34, expression data using empirical Bayes methods. Biostatistics 8, 118–127 (2007). 726–737 (2016). 53. Rosenbloom, K. R. et al. ENCODE Data in the UCSC Genome Browser: year 5 update. Nucleic Acids Res. 41, D56–D63 (2013). Acknowledgements 54. Davis, C. A. et al. The Encyclopedia of DNA elements (ENCODE): Data portal This work was supported by the Post-Genome Technology Development Program update. Nucleic Acids Res. 46, D794–D801 (2018). [10067758, Business model development driven by clinico-genomic database for preci- 55. Liberzon, A. et al. The molecular signatures database hallmark gene set sion immuno-oncology] funded by the Ministry of Trade, Industry and Energy (MOTIE, collection. Cell Syst. 1, 417–425 (2015). Korea). H.J. was supported by a grant from the National Research Foundation of Korea 56. Thorsson, V. et al. The immune landscape of cancer. Immunity 48, 812–830 funded by the Korean Government (NRF-2018R1A6A3A01010889). (2018). e14. 57. Charoentong, P. et al. Pan-cancer immunogenomic analyses reveal genotype- Author contributions immunophenotype relationships and predictors of response to checkpoint H.J. performed all data analyses and wrote the manuscript. H.S.K. generated all data and blockade. Cell Rep. 18, 248–262 (2017). helped with manuscript writing. J.Y.K. participated in data analyses. J.-M.S., J.S.A., 58. Evrony, G. D. et al. Cell lineage analysis in human brain using endogenous M.-J.A., K.P. and M.E. helped with cohort data generation and analyses. S.-H.L. and retroelements. Neuron 85, 49–60 (2015). J.K.C. conceived and supervised the study. 59. Mayer, J. et al. A revised nomenclature for transcribed human endogenous retroviral loci. Mob. DNA 2, 7 (2011). 60. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows- Additional information Wheeler transform. Bioinformatics 25, 1754–1760 (2009). Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 61. Chen, R., Im, H. & Snyder, M. Whole-exome enrichment with the agilent 019-12159-9. sureselect human all exon platform. Cold Spring Harb. Protoc. 2015, 626–633 (2015). Competing interests: The authors declare no competing interests. 62. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, Reprints and permission information is available online at http://npg.nature.com/ 1297–1303 (2010). reprintsandpermissions/ 63. Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 1–4 (2018). Peer review information Nature Communications thanks the anonymous reviewer(s) for 64. Sherry, S. T. et al. dbSNP: the NCBI database of genetic variation. Nucleic their contribution to the peer review of this work. Peer reviewer reports are available. Acids Res. 29, 308–311 (2001). 65. Yang, H. & Wang, K. Genomic variant annotation and prioritization with Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in ANNOVAR and wANNOVAR. Nat. Protoc. 10, 1556 (2015). published maps and institutional affiliations. 66. Talevich, E., Shain, A. H., Botton, T. & Bastian, B. C. CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing. PLOS Comput. Biol. 12, e1004873 (2016). Open Access This article is licensed under a Creative Commons 67. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, Attribution 4.0 International License, which permits use, sharing, (15–21 (2013). adaptation, distribution and reproduction in any medium or format, as long as you give 68. Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq appropriate credit to the original author(s) and the source, provide a link to the Creative data with or without a reference genome. BMC Bioinf. 12, 323 (2011). Commons license, and indicate if changes were made. The images or other third party 69. Aryee, M. J. et al. Minfi: A flexible and comprehensive Bioconductor package material in this article are included in the article’s Creative Commons license, unless for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, indicated otherwise in a credit line to the material. If material is not included in the 1363–1369 (2014). article’s Creative Commons license and your intended use is not permitted by statutory 70. Tian, Y. et al. ChAMP: updated methylation analysis pipeline for Illumina regulation or exceeds the permitted use, you will need to obtain permission directly from BeadChips. Bioinformatics 33, (3982–3984 (2017). 71. Chakravarthy, A. et al. Pan-cancer deconvolution of tumour composition the copyright holder. To view a copy of this license, visit http://creativecommons.org/ using DNA methylation. Nat. Commun. 9, 3220 (2018). licenses/by/4.0/. 72. Yoshihara, K. et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612(2013). © The Author(s) 2019 12 NATURE COMMUNICATIONS | (2019)10:4278 | https://doi.org/10.1038/s41467-019-12159-9 | www.nature.com/naturecommunications