Abstract
Background:
Epigenetic changes of lung adenocarcinoma (LUAD) have been reported to be a relevant factor in tumorigenesis and cancer progression. However, the molecular mechanisms responsible for DNA methylation patterns in the tumor immune-infiltrating microenvironment and in cancer immunotherapy remain unclear.
Methods:
We conducted a global analysis of the DNA methylation modification pattern (DMP) and immune cell-infiltrating characteristics of LUAD patients based on 21 DNA methylation regulators. A DNA methylation score (DMS) system was constructed to quantify the DMP model in each patient and estimate their potential benefit from immunotherapy.
Results:
Two DNA methylation modification patterns able to distinctly characterize the immune microenvironment characterization were identified among 513 LUAD samples. A lower DMS, characterized by increased CTLA-4/PD-1/L1 gene expression, greater methylation modifications, and tumor mutation burden, characterized a noninflamed phenotype with worse survival. A higher DMS, characterized by decreased methylation modification, a greater stromal-relevant response, and immune hyperactivation, characterized an inflamed phenotype with better prognosis. Moreover, a lower DMS indicated an increased mutation load and exhibited a poor immunotherapeutic response in the anti-CTLA-4/PD-1/PD-L1 cohort.
Conclusion:
Evaluating the DNA methylation modification pattern of LUAD patients could enhance our understanding of the features of tumor microenvironment characterization and may promote more favorable immunotherapy strategies.
Introduction
The global DNA methylation is strongly associated with growth selection and uncontrolled cell proliferation in multiple cancer types (1, 2). Among the epigenetic mechanisms of the mammalian genome, DNA methylation is catalyzed by series of DNA methyltransferases catalyzing the transfer of a methyl group from S-adenyl methionine to the C5 position of a cytosine residue to form 5-methylcytosine at CpG sites (3). The expression and function of those methyltransferases has been reported to be partly responsible for immunomodulation and might have an impact on DNA methylation modification patterns in human cancers (4, 5).
Lung adenocarcinoma (LUAD) is the most common form of lung cancer and is responsible for the cancer deaths worldwide (6). However, the overall survival time of LUAD patients remains short despite the enormous research efforts for the development of effective diagnostic techniques and therapeutics. Recently, as a result of an increasing number of studies focusing on tumor immune cell infiltrations and DNA epigenetic modification, crucial immune cell types and methyltransferase subsets in tumor growth, metastasis, and outcome have gradually been identified (7–10). The integrated analysis of the composition of the immune cells in both LUAD tumors and paired normal adjacent tissues has revealed that a deeper exploration of immune signatures or biomarkers in the tumor microenvironment (TME) could play an essential role in revealing the potential oncogenic mechanism in LUAD (11). Moreover, the immunotherapeutic efficacy of several immune checkpoint inhibitors (ICIs) (anti-CTLA4/PD-L1/PD-1, 8) has widely assessed and has achieved a notable response in tumor treatment, including LUAD (12–14). Further studies have also proposed that the identification of altered epigenetic methylation patterns may represent a valuable diagnostic approach toward novel therapeutic strategies for preventing and treating LUAD. For instance, Yang et al. reported that the downregulation of methyltransferases DNMT3A and MBD4 could promote ALDH2 expression, reducing the probability of bone metastasis in LUAD patients (15). Forloni et al. found that the inhibition of the demethylase TET oncogene family member 1 (TET1) could induce epigenetic silencing of antitumor genes in the oncogenesis EGFR signaling pathway, indicating that dysregulated DNA methylation probably played a major role in tumorigenesis (16). However, the global profiles of DNA methylation regulators based on the correlation between immune microenvironment and immunotherapy of LUAD samples have not been fully evaluated (12, 17).
In this study, we collected and integrated the clinical information and genomic data of 513 patients from The Cancer Genome Atlas (TCGA) LUAD cohort and comprehensively evaluated the TME characteristics represented by distinct patterns of DNA methylation modifications. We finally identified two independent DNA methylation patterns (DMPs) by unsupervised consensus clustering the expression of 21 DNA methylation regulator-related genes, which we defined as methylation regulators. The immune cell-infiltrating properties of two DMPs were highly consistent with immune-noninflamed or immune-inflamed phenotype, respectively. We then constructed a DNA methylation score (DMS) system to investigate the efficacy of DNA methylation modification patterns in individual patients and estimated their immunotherapeutic value in several clinical trials.
Materials and Methods
Data Collection and Preparation
The workflow of this study is presented in Figure S1A. Gene expression data and complete clinical annotations of LUAD samples were retrospectively collected from publicly available datasets of the Gene-Expression Omnibus (GEO) and TCGA databases. A total of six GEO lung adenocarcinoma cohort somatic mutation, copy number variation (CNV), and clinical data, including tumor stage, age, sex, and overall survival times/states were obtained from TCGA databases (GSE116959, GSE58772, GSE99995, GSE68571, GSE68465, and GSE26939) and one TCGA-LUAD cohort were enrolled for further analysis. For the GEO microarray data, the normalized matrix files were directly downloaded. As for TCGA cohort, the RNA sequencing data (FPKM value) for gene expression was downloaded from the Genomic Data Commons and then transformed into transcripts per kilobase million (TPM) format. Batch effects among GEO datasets were corrected using the “ComBat” algorithm of the sva package. The baseline information of all eligible datasets is summarized in Table S1.
Consensus Molecular Clustering of 21 DNA Methylation Regulators
We collected DNA methylation regulator-related studies, and a total of 21 regulator genes (including “writers’: DNMT1, DNMT3A, DNMT3B; “easers”: TET1, TET2, TET3; and “readers”: MBD1, MBD2, MBD3, MBD4, TDG, SMUG1, UHRF1, UHRF2, ZBTB4, ZBTB24, ZBTB33, ZBTB38, NSUN2, MGMT, DMAP1) were extracted for DNA methylation modification pattern identification (Table S2). The ConsensusClusterPlus algorithm was employed to conduct unsupervised clustering of individual tumor samples with the gene expression profiles of 21 DNA methylation regulators (18). The cluster assignments were stable when k = 2.
Gene-Set Variation Analysis and Functional Annotation
We performed Gene-Set Variation Analysis (GSVA) enrichment analysis with “GSVA” R packages to study the differed biological pathway between different DNA methylation modification patterns in cancer samples (19). The Hallmarker gene-set of “c2.cp.kegg.v7.4.symbols” was retrieved from the MSigDB database (20). Adjusted p-values <0.05 were considered statistically significant for GSVA analysis. Gene ontology (GO) and pathway annotation for DNA methylation pattern-related genes were performed using the R package “clusterProfiler” with a cutoff value of p < 0.05.
Immune Cell Infiltration Estimation
The Single Sample Gene-Set Enrichment Analysis (ssGSEA) was implemented to determine variations in immune leukocyte subtype abundance between different DMP clusters using the R package “GSEAbase.” Subsequently, the abundances of 22 immune cell types for each tumor specimen were further identified by estimating relative subsets of RNA transcripts (CIBERSORT; https://cibersort.stanford.edu/) using the gene expression profile of LUAD cancer (21).
Identification of Differentially Expressed Genes and DMS Construction
We used three R packages (“limma,” “edgeR,” and “Deseq2”) to identify the DNA methylation modification-related differentially expressed genes (DEGs) across distinct DMP phenotypes. Univariate Cox model analysis was performed to calculate their association with overall survival and to extract prognostic DEGs to construct a scoring system. We then conducted principal component analysis (PCA) using the identified DEG prognostic genes, and PCA 1–2 components were selected to act as signature scores to construct DNA methylation-relevant gene signature, which we defined as the DMS. This method mainly focused on evaluating the score for each patient in the dataset with the largest group of well-correlated (or anticorrelated) genes. We constructed the DMS using an algorithm similar to the GGI: DMS = Σ(PC1i + PC2i), where i stood for the relevant gene expression of the selected set (22, 23).
Immunotherapy Dataset Collection
We searched the immunotherapeutic characteristics of the TCGA-LUAD patients using the Cancer Immunome Database (TCIA), which is a publicly available dataset containing corresponding clinical pathology information (24). The immunotherapeutic GEO datasets were included in this study: GSE135222 (anti-PD-1/PD-L1 treatment) (25) and GSE91061 (anti-CTLA-4/PD-1 treatment) (26). We also retrieved IMvigor210 datasets containing data on atezolizumab treatment and extracted the relative gene expression profiles and clinical notes using the R “IMvigor210CoreBiologies” package (27). The raw count data were then transformed into TPM value format.
Results
Landscape of DNA Methylation Regulators in LUAD
In this study, 21 DNA methylation regulators were collected from the published data. Figure 1A displays the dynamic reversible process of RNA methylation mediated by regulators and their corresponding potential biological functions. We further analyzed the genetic alterations of DNA methylation regulators in LUAD. Among the 561 cases, 123 samples (21.93%) harbored genetic alterations of DNA methylation regulators, primarily including nonsense or missense mutations. TET1 (4%) had the highest mutation frequency, followed by DNMT3A, TET3, and DNMT3B (Figure 1B). A mutation cooccurrence pattern across all DNA methylation regulators was examined, and significant gene mutation relationships were identified between DNMT1 and ZBTB4, MBD4 and NSUN2, and DNMT1 and UHRF1 (Figure S1B). Further analysis of CNV alteration frequency in the 21 regulators showed a prevalence of CNV mutations. Among these, NSUN2, DMAP1, SMUG1, DNMT3B, ZBTB33/38, and DNMT3A showed relatively higher amplification frequency in terms of CNV, while MBD1/2/3, ZBTB4/24, UHRF1, TDG, and TET2 had a widespread frequency of copy number deletion variants (Figure 1C). The location of CNV among the DNA methylation regulators on chromosomes is shown in Figure 1D. To investigate whether these genetic variants influenced the gene expression of regulators in cancer samples, we compared the expression of the 21 DNA methylation regulators between normal and LUAD samples and found that regulators with significant genetic mutation as well as CNV amplification were significantly higher expressed in cancer groups (Figures 1B–E), indicating that the regulator gene mutation and the alternation of CNV could be risk factors leading to the multiperturbation of the translation of DNA methylation regulators. Univariable and multivariable Cox proportional hazards regression analyses were performed to evaluate the association between regulators and patient survival. The Forest plot and Kaplan-Meier curve showed that DMAP1 and MGMT were upregulated in tumor specimens and could be considered protective factors associated with prolonged survival time, while ZBTB38 and MBD2/3 downregulation in the tumor group was recognized as a risk factor with worse overall survival (Figures S1C–H), indicating that the unbalance expression of DNA methylation regulators could contribute to tumorigenesis. These analyses demonstrated that the crosstalk between DNA methylation regulators in the genomic and transcriptomic landscape has a potential role in LUAD occurrence and complexity.
Figure 1
DNA Methylation Modification Patterns in TCGA Cohort
To determine the pattern of DNA methylation modification mediated by 21 regulators, the TCGA-LUAD cohort and their corresponding clinical data were prepared for analysis (Table S3). A comprehensive landscape of the interactions and the prognostic significance of the 21 regulators was visualized using the conetwork plot (Figure 2A). We dissected the relationship among those regulators and found that most regulators showed a positive correlation with each other, which is consistent with the above analysis. We also demonstrated that the expression of UHRF1, SMUG1, TDG, and MBD1/2/3 presented the most significant positive correlation with other regulators and might be the risk factors for carcinoma formation. Notably, the readers, such as MGMT and ZBTB4, showed a remarkably negative relationship with the other regulators. Thus, these findings indicated that crosstalk among DNA methylation regulators may play a critical role in the onset of LUAD.
Figure 2
We continued to conduct unsupervised clustering based on the expression of 21 DNA methylation regulators to identify DNA methylation modification patterns in LUAD samples, and two distinct clusters were accordingly obtained via the ConsensusClusterPlus package in R software, including 246 cases in cluster A and 267 cases in cluster B (Figure 2B). We renamed those clusters as DNA methylation pattern (DMP-A and DMP-B) (Table S4). Survival analysis for the two DNA methylation modification clusters showed that DMP-B had a relatively better prognosis in the TCGA cohort (Figure 2C). In addition, we conducted heatmap analysis of the relative expression of 21 DNA methylation regulators between distinct DNA methylation modification patterns and observed that most regulators were markedly elevated in the DMP-A cluster (Figure 2B). Principal component analysis (PCA) also indicated that two DMP clusters were completely segregated, indicating that they could be easily distinguished via the gene expression pattern of the 21 regulators (Figure 2D). The heatmap analysis of the tumor methylation status determined that the DMP-B cluster presented a low-methylation epigenotype with worser clinical outcomes (Figure 2F). These results demonstrated that altered expression of DNA methylation regulators could affect patient survival by regulating tumor methylation levels, which contributes to the high heterogeneity of LUAD.
TME Cell Infiltration Characteristics of the Two Distinct DMPs in the TCGA Cohort
To further explore differences in biological behaviors among the two DNA methylation modification patterns, we performed GSVA enrichment analysis. As shown in Figure 2E and Table S5, DMP-A was mainly presented a CD4+ immune regulator and stromal activation phenotype, associated with many related pathways, such as the cell cycle, DNA replication, Toll-like receptor, and Nod-like receptor signaling pathways in cancer. Whereas, DMP-B represented enriched pathways associated with immune/inflammation hyperactivation including the activation of the PPAR signaling pathway and the autoimmune response, including asthma, autoimmune thyroid disease, allograft rejection, and graft-versus-host disease. Furthermore, we conducted a TME cell infiltration analysis using the ssGSVA algorithm and demonstrated that DMP-A showed a relatively lower proportion of immune cell infiltration and exhibited a shorter survival time (Figures 2B, 3A). However, DMP-B did present a matching survival advantage and significantly elevated innate immune cells including activated B cells, MDSCs, macrophages, mast cells, natural killer cells, and activated CD8+ T cells and so on. The CIBERSORT analysis also indicated that DMP-B was markedly enriched in immunoactive cells, such as M2 macrophages, memory B cells, plasma cells, and Tregs, which is coherently indicative of a better prognosis (Figures 2B, 3B). Previous studies have reported that the immune-excluded phenotype in tumors also presented immune cells limited to the stroma surrounding the tumor cell nests (28). Therefore, we suspected that the stromal activation and low immune cell density in DMP-A pattern caused the suppression of antitumor response, which was also validated in the GEO LUAD cohort (GSE116959, GSE58772, and GSE99995) (Figure S2I). Altogether, these results indicated that LUAD patients could be classified into two immune phenotype groups, among which the DMP-A pattern represented a noninflamed phenotype (immune-excluded or immune-deserted trait) characterized by stromal activation and weakened immune regulator status, whereas the DMP-B pattern represented an immune-inflamed phenotype characterized by immune/inflammation hyperactivation and multicomponent immune cell infiltration in the TME.
Figure 3
DMP Phenotype-Related DEGs and Functional Annotation
To further explore the characteristics of the DMP in the LUAD cohort, a total of 2,720 DNA methylation regulator pattern related to the DEGs were obtained and subjected to conduct enrichment analysis. The KEGG pathway related to the cell cycle and the p53 signaling pathway were also significantly enriched in DMP-related gene set (Figure 3C). The results of GO terms showed that those genes were mainly enriched in the cell cycle process and DNA replication (Figure 3D). We subsequently selected 792 prognostic genes (p < 0.05) via univariate Cox model analysis (Table S6). Based on the expression of the 792 representative genes, unsupervised clustering analysis was performed to classify LUAD patients into three distinct clusters, namely, DNA methylation gene cluster I/II/III (Figures 4A, B; Figures S2A–E). The heatmap and survival analysis showed that patients in gene cluster III were characterized by earlier clinical stage and better prognosis, while patients in gene cluster II were characterized by advanced staging and worser survival outcome (Figures 4B, C). In addition, multivariate Cox regression analysis indicated that the overall survival between the gene cluster gene signatures remained statistically significant after considering multiple factors, including age, tumor stage (T and N), and sex (Figure 4D), indicating that the identified gene clusters could represent independent factors of prognosis. We also noted that the gene expression of the 21 regulators exhibited significant differences among the three gene clusters, and in particularly patients in gene cluster III were related to the higher expression of protective factors (TET2/MGMT), which is in line with the findings regarding DNA methylation modification patterns (Figure 4E). The patients in gene cluster III also presented higher immune cell infiltration, indicating a phototype of immune activation. While the levels of almost all immune cells in gene cluster II were relatively lower than in other groups, these presented a phenotype of immune suppression (Figure 4F). We further evaluated the DNA methylation status of the gene clusters, and results indicated that the gene cluster II exhibited a hypermethylation status, whereas the gene cluster III presented a weaker methylation status (Figure 4G). All these results also confirmed that the perturbed DNA methylation regulator genetic expression was highly correlated with different immune responses, genetic methylation levels, and clinical outcomes and contributed to LUAD heterogeneity, which was consistent with the previous findings.
Figure 4
Construction of DMS and Exploration of Clinical Relevance
To accurately evaluate the DNA methylation status of individual patients with LUAD, we developed a scoring scheme termed the DMS based on the expression of 792 DNA methylation-related genes.
The DMP-B pattern had a higher level of DMS compared with patients in DMP-A (Figure 5A). Notably, gene cluster III showed the highest level of DMS, followed by gene clusters I and II (Figure 5B). With an optimum cutoff value of 2.17 determined by the survminer package, we divided LUAD patients into high- and low-DMS groups. The alluvial diagram summarized the attribute changes of patients according to the DNA methylation regulator pattern, gene cluster, and DMS groups (Figure 5C). We then examined the correlation between biological processes, immune cell infiltration, and the level of DMS signatures using Spearman’s analysis. The DMS signatures were significantly negatively correlated with cell cycle and DNA/RNA repair signatures but were positively correlated with immune activation and stromal-relevant pathways (Figure 5D; Figure S2F). In addition, we calculated the distribution differences of tumor genome somatic mutations between the DMS groups using the R “maftools” package. As we expected, the low-DMS group with shorter survival time presented significantly higher tumor mutational profiles than the high-DMS group, consistent with recent findings that gene alterations correlated highly with tumor invasion and cell proliferation. Furthermore, survival analysis revealed that the DMS-low patients were significantly correlated with a worse prognosis in the TCGA cohort (p < 0.001, Figure 5E), which was further validated by the ROC curves (AUC = 0.652, Figure 5F). The multivariate Cox regression model confirmed the DMS factor could stand as an independent prognostic biomarker for evaluating patient outcomes in the TCGA-LUAD cohort (HR, 2.03 [95% confidence interval, 1.34–3.08], p < 0.001, Figure 5G). To further verify this DMS model, we also performed multivariate Cox regression model and prognosis analysis using meta-GEO cohorts (GSE68465, GSE68571, and GSE26939), which provided further support that DMS was a significant prognostic factor for predicting patient outcomes (Figures S2G, H). The analysis of tumor mutation burden (TMB) confirmed that the low-DMS cluster was significantly correlated with a higher TMB (Figures 5H, I). The mutational landscape of significantly mutated genes showed that most genes (including TP53 and TTN) had higher somatic alteration rates in the lower-DMS group, whereas KRAS had a higher somatic alteration rate in the higher-DMS group. We further observed that there was a markedly negative correlation between the TMB score and DMS (R = −0.48, p < 0.001), demonstrating the crosstalk between the DMS and genetic mutation evaluation (Figure 5J). Meanwhile, the methylation heatmap analysis revealed that the higher-DMS group exhibited hypomethylation compared with the lower-DMS group (Figure S2J). These results demonstrated that the DMS could represent the DNA methylation modification patterns and comprehensively could reflect genomic variation modifications, as well as effectively predict the prognosis of LUAD patients.
Figure 5
The Role of DMS in Predicting Immunotherapy
Significant progress has been made to identify the effective immune-related signatures that correlated with response to antitumor therapy, particularly represented by TMB and ICIs. Our analysis found that the gene expression of the CTLA-4/PD-1/PD-L1 in TCGA-LUAD cohort was significantly increased in the DMS-low group (Figures 6A–C). Considering the high correlation across DMS levels with the immune response, we explored whether the DMS system could predict the patients’ response to ICI treatment in two independent immunotherapy cohorts treated with CTLA-4/PD-1/PD-L1 antibody inhibitors. In the TCIA-LUAD cohort (anti-CTLA-4/PD-1 treatment), patients classified as the DMS-high group exhibited a significantly clinical response (Figures 6D–G), indicating the immunotherapeutic benefits of CTLA-4/PD-1 antibody treatment in DMS-high patients. In the anti-CTLA-4/PD-1/PD-L1 GEO cohort (GSE91061 and GSE135222), higher DMS in tumor patients was associated with a stronger immune response and higher therapeutic benefits (Figures 6H, I). The progression-free survival analysis showed that higher-DMS patients exhibited prolonged survival (Figure 6I). A similar result was also identified in the IMvigor210 cohort (anti-PD-L1 treatment) (Figures 6J–M). Those results demonstrated that the DMS signature system was strongly associated with the tumor immune response and might contribute to predicting the efficacy value of the anti-CTLA-4/PD-1/PD-L1 immunotherapy.
Figure 6
Discussion
Increasing evidence has indicated that aberrant DNA methylation might increase genome instability by silencing of notable antioncogenes by methylated modifications, resulting in TME alterations, CNV, and biological process conversion (29, 30). However, the global modulation of DNA methylation modification in the immune contexture of LUAD patients remains to be comprehensively recognized. Thus, the identification of distinct DNA methylation modification patterns in the tumor immune microenvironment could provide a useful evaluation of the correlation between DNA methylation on the immune-related response and may assist in achieving more effective immunotherapy resolution.
Herein, following unsupervised clustering analysis of the gene expression of 21 regulators, we defined two different DMPs characterized by distinct immune phenotypes in the tumor immune microenvironment. DMP-A–clustered patients presented with hypermethylation had a lower proportion of immune cells, characterized by the suppression of the immune cell-infiltrating response in tumor cells and a noninflamed (immune-excluded/deserted) phenotype, which corresponded to a worser survival prognosis. In contrast, the DMP-B–clustered individuals presented with hypomethylation was characterized by the activation of abundant immune effector/pathways and the presence of multiple immune cell infiltrations in the TME, which exhibited an immune-inflamed phenomenon, corresponding to significant survival benefits. Analysis of the tumor immunophenotype in patients with solid tumors has also widely supported the view that the immune environment plays a central role in the pathogenesis of tumorigenesis and metastasis and harbored important clinical implications for effective immunotherapy outcomes (27, 31–33). The available literature to date has reported that the immune-excluded phenomenon is largely reflected by hyper-invigorating stromal cells surrounding tumor cell nests, categorized as an immune-inflamed phenotype (34, 35). Stroma with lower expression of immune markers (tumor-infiltrating lymphocytes and CD8+ T cells) could suppress the antitumor signatures and interfere with the penetration of immune cell infiltration into the tumor parenchyma, consistent with the poorer survival in the DMP-A immune phenotype (36–38). The immune-inflamed phenotype was represented by a higher density of immune T cells, enrichment of abundant immune signal pathways, and the presence of preexisting immune infiltration in an antitumor microenvironment (34, 39). Recapitulating our analysis of the TME immune-infiltrating characteristics using the ssGSEA and CIBERSORT algorithm with our proposed DMP clusters confirmed the above definition of immune classification. In addition, previous reports have also shown that DNA hypomethylation is associated with immune signaling activation (40), which further supported our findings in this study.
To further elucidate the characteristics of transcriptome traits in DMP, the DEGs of distinct DNA methylation modification patterns were obtained and selected as the DNA methylation signature genes, which were concurrently found to be significantly associated with the cell cycle and DNA replication biological behaviors. Based on the mRNA expression levels of these DEGs, three DNA methylation gene clusters with different survival prognosis and immune-infiltrating features were identified, indicating there were indeed different clinical immune phenotypes associated with LUAD (11, 41). In addition, we established a methylation-based scoring system to quantify the DNA methylation modifications of individual patients, which helped to define the DMP signature, thus yielding greater scientific insight into the complex mechanisms of tumorigenesis and progression (27, 31, 32). As expected, the DMP-A group was associated with gene cluster I/II and was characterized by the presence of immune suppression and a higher methylation phenotype and a lower DMS, which corresponded to a shorter survival time. The multivariate analysis of the LUAD cohort confirmed the DMS system could be an independent factor for patient prognosis. Altogether, the above results indicated that there was indeed a higher heterogeneity of DNA methylation modifications in the tumor immune/alteration microenvironment of LUAD (42).
Consistent with the previous analysis, the DMS, which was negatively correlated with TMB, was found to be markedly enriched in the immune activation and stromal-relevant pathways, underlining the pivotal role of stromal immune activation in resistance to immune checkpoint therapeutics (43, 44). Our study revealed that cancer drug resistance was not only accompanied by the hyperactivation of various stromal pathway associated with TGF beta, MAPK, and Wnt signaling pathway but also correlated to gene repair and cell cycle processes. For instance, the inhibition of the TGF beta pathway was reported to induce durable responses to PD-1/PD-L1 blockade in tumor models (45, 46). Previous studies also demonstrated that the cancer epithelial-mesenchymal transition is frequently activated by TGF-β, Wnt, Notch, and MAPK signaling pathways, which are the major factors promoting metastasis and notorious invasion of cancer cells (47, 48). Furthermore, recent studies searching for the effective immunotherapies in cancer and found that ICIs were successful cancer treatments, particularly in metastatic urothelial cancer and melanoma where anti-CTLA-4/PD-1/PD-L1 antibodies have found widespread application (27, 49, 50). However, using our data, we evaluated the therapeutic value of the DMS in the immune checkpoint (CTLA-4/PD-1/PD-L1) treatment cohorts and showed the opposite clinical benefit to immune inhibitors. Compared with the lower-DMS group, the higher-DMS group, which presented a smaller proportion of genetic mutations and methylation, exhibited a stronger antitumor immune response and benefit of ICI treatment. Meanwhile, TP53 and TTN mutations in the lower-DMS subgroup showed a larger proportion of mutation rates compared with the higher-DMS subgroup, whereas the KRAS mutation rate increased in the high-DMS subgroup. TP53 and KRAS are prevalent oncogenic drivers in most tumor types, and their cooccurring mutations result in the upregulation of tumor immunogenicity and immune tolerance/escape in response to PD-1 blockade immunotherapy in LUAD (51). Furthermore, TTN and TP53 mutations may exhibit a synergistically prognostic benefit in various lung cancers except LUAD (52). Those findings supported that the notion that the DMS signature of DNA methylation regulator patterns, combined with gene mutation signals, as promising to predict the efficacy of anti-CTLA-4/PD-1/PD-L1 immune checkpoint blockade therapy, and could contribute to guide more effective strategies for precision immunotherapy in cancer individuals. Therefore, the DMS signature could stand as an important biomarker to estimate the benefit of antitumor immune response to cancer treatment (53, 54).
However, there were some limitations in this study. Firstly, we used only 21 DNA methylase-related regulators to construct the model system; newer methylation regulators should be included in further analyses. Second, additional clinical data are needed to validate the efficacy of the immunotherapy and to evaluate the accuracy of the DMP and DMS system.
Conclusion
We identified two DNA methylation modification patterns in LUAD patients based on 21 methylase-related regulators and constructed a methylation profile score model for individual patients. The TME disparity between the distinct DMPs highlighted the potential complexity and heterogeneity of LUAD formation, casting light on multiple processes fueling tumor evolution, immune infiltration, and drug resistance. The systematic analysis of the immune and/or clinical characteristics in the DMP and DMS will contribute to enhance a deeper understanding of the tumor immune-infiltrating microenvironment and will promote the development of effectively targeted immunotherapies.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Statements
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.
Author contributions
DY designed the experiment, analyzed the data, drafted the manuscript, and coordinated the study. ZW revised the paper. YW, FC, LY, and SZ performed the table and figure processing. RT and JL supervised the study. All authors contributed to the article and approved the submitted version.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.734873/full#supplementary-material
Supplementary Figure 1Correlation and prognostic value of 21 DNA methylation regulators and overview of study design. (A) Overview of this study flow. (B) The cooccurrence and mutual exclusion of 21 regulators. Cooccurrence, blackish green; mutual exclusive, yellowish brown. (C, D) Clinical prognostic significance of expression of DNA regulators by univariate (C) and multivariate (D) Cox regression analyses. (E–H) Survival analysis of the regulators in TCGA cohort, including MBD3, DMAP1, and ZBTB38.
Supplementary Figure 2Unsupervised clustering of 586 DNA methylation gene signatures and construction of the DMS system. (A–C) The consensus matrix (CM) plot of DNA methylation gene signature for k = 3 to k = 5. (D) The distribution of the consensus index cumulative distribution function (CDF) curves in the LUAD cohort for each k. (E) Tracking plot of subclusters at each k. (F) Correlation between DMS and immune-infiltrating trait in TCGA samples by the Spearman’s analysis. (G) Subcluster analysis determining the clinical prognostic ability of DMS in the meta-GEO cohort by multivariate Cox regression. The length of the horizontal line represented the 95% confidence interval for each group. (H) Survival analysis of patient between low- or high-DMS subclusters in the meta-GEO cohort. (I) The abundance of immune cell infiltration in different DMP groups using the ssGSVA algorithm in the meta-GEO cohort. (J) The methylation heatmap analysis of low- or high-DMS group.
Supplementary Table 1Clinical information of LUAD cohorts from GEO/TCGA in this study.
Supplementary Table 2Summary of DNA methylation modification regulators.
Supplementary Table 3Clinical annotation and DNA methylation regulators of individual patient in TCGA-LUAD cohort.
Supplementary Table 4Information of samples clustering in TCGA-LUAD cohorts.
Supplementary Table 5Enrichment score of KEGG pathways in TCGA-LUAD cohorts.
Supplementary Table 6Prognostic DEGs related to DNA methylation modification patterns using the univariate Cox regression model.
Abbreviations
CNV, copy number variation; DEGs, differentially expressed genes; DMP, DNA methylation modification pattern; DMPs, DNA methylation patterns; DMS, DNA methylation score; GEO, Gene-Expression Omnibus; GO, Gene Ontology; GSVA, gene-set variation analysis; ICIs, immune checkpoint inhibitors; LUAD, lung adenocarcinoma; PCA, principal component analysis; ssGSEA, single sample gene-set enrichment analysis; TCGA, The Cancer Genome Atlas; TCIA, Cancer Immunome Database; TME, tumor microenvironment; TMB: tumor mutation burden; TPM, transcripts per kilobase million.
References
1
RahmaniMTalebiMHaghMFFeiziASolaliS. Aberrant DNA Methylation of Key Genes and Acute Lymphoblastic Leukemia. Ritorno Al Numero (2018) 97:1493–500. doi: 10.1016/j.biopha.2017.11.033
2
LasseigneBNBrooksJD. The Role of DNA Methylation in Renal Cell Carcinoma. Mol Diagn Ther (2018) 22(4):431–42. doi: 10.1007/s40291-018-0337-9
3
MooreLDLeTFanG. DNA Methylation and its Basic Function. Neuropsychopharmacology (2012) 38(1):23–38. doi: 10.1038/npp.2012.112
4
LeeSKimHSRohKHLeeBCShinTHYooJMet al. DNA Methyltransferase Inhibition Accelerates the Immunomodulation and Migration of Human Mesenchymal Stem Cells. Sci Rep (2015) 5:8020. doi: 10.1038/srep08020
5
SigalottiLCovreAFrattaEParisiGColizziFRizzoAet al. Epigenetics of Human Cutaneous Melanoma: Setting the Stage for New Therapeutic Strategies. J Transl Med (2010) 8:56. doi: 10.1186/1479-5876-8-56
6
BrayFFerlayJSoerjomataramISiegelRLTorreLAJemalA. Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492
7
BruniDAngellHKGalonJ. The Immune Contexture and Immunoscore in Cancer Prognosis and Therapeutic Efficacy. Nat Rev Cancer (2020) 20(11):662–80. doi: 10.1038/s41568-020-0285-7
8
HashimotoHVertinoPMChengX. Molecular Coupling of DNA Methylation and Histone Methylation. Epigenomics (2010) 2(5):657–69. doi: 10.2217/epi.10.44
9
RasmussenKDHelinK. Role of TET Enzymes in DNA Methylation, Development, and Cancer. Genes Dev (2016) 30(7):733–50. doi: 10.1101/gad.276568.115
10
ThompsonJJKaurRSosaCPLeeJHKashiwagiKZhouDet al. ZBTB24 is a Transcriptional Regulator That Coordinates With DNMT3B to Control DNA Methylation. Nucleic Acids Res (2018) 46(19):10034–51. doi: 10.1093/nar/gky682
11
GilletteMASatpathySCaoSDhanasekaranSMVasaikarSVKrugKet al. Proteogenomic Characterization Reveals Therapeutic Vulnerabilities in Lung Adenocarcinoma. Cell (2020) 182(1):200–25.e235. doi: 10.1016/j.cell.2020.06.013
12
ZhangLZhangZYuZ. Identification of a Novel Glycolysis-Related Gene Signature for Predicting Metastasis and Survival in Patients With Lung Adenocarcinoma. J Transl Med (2019) 17(1):423. doi: 10.1186/s12967-019-02173-2
13
DubeyDDavidWSAmatoAAReynoldsKLClementNFChuteDFet al. Varied Phenotypes and Management of Immune Checkpoint Inhibitor-Associated Neuropathies. Neurology (2019) 93(11):e1093–103. doi: 10.1212/wnl.0000000000008091
14
FarukiHMayhewGMSerodyJSHayesDNPerouCMLai-GoldmanM. Lung Adenocarcinoma and Squamous Cell Carcinoma Gene Expression Subtypes Demonstrate Significant Differences in Tumor Immune Landscape. J Thorac Oncol (2017) 12(6):943–53. doi: 10.1016/j.jtho.2017.03.010
15
YangMWangALiCSunJYiGChengHet al. Methylation-Induced Silencing of ALDH2 Facilitates Lung Adenocarcinoma Bone Metastasis by Activating the MAPK Pathway. Front Oncol (2020) 10:1141. doi: 10.3389/fonc.2020.01141
16
ForloniMGuptaRNagarajanASunLSDongYPirazzoliVet al. Oncogenic EGFR Represses the TET1 DNA Demethylase to Induce Silencing of Tumor Suppressors in Cancer Cells. Cell Rep (2016) 16(2):457–71. doi: 10.1016/j.celrep.2016.05.087
17
XuFHeLZhanXChenJXuHHuangXet al. DNA Methylation-Based Lung Adenocarcinoma Subtypes can Predict Prognosis, Recurrence, and Immunotherapeutic Implications. Aging (Albany NY) (2020) 12(24):25275–93. doi: 10.18632/aging.104129
18
WilkersonMDHayesDN. ConsensusClusterPlus: A Class Discovery Tool With Confidence Assessments and Item Tracking. Bioinformatics (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170
19
HänzelmannSCasteloRGuinneyJ. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7
20
LiberzonABirgerCThorvaldsdóttirHGhandiMMesirovJPTamayoP. The Molecular Signatures Database (MSigDB) Hallmark Gene Set Collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004
21
NewmanAMLiuCLGreenMRGentlesAJFengWXuYet al. Robust Enumeration of Cell Subsets From Tissue Expression Profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
22
ZengDLiMZhouRZhangJSunHShiMet al. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res (2019) 7(5):737–50. doi: 10.1158/2326-6066.Cir-18-0436
23
SotiriouCWirapatiPLoiSHarrisAFoxSSmedsJet al. Gene Expression Profiling in Breast Cancer: Understanding the Molecular Basis of Histologic Grade to Improve Prognosis. J Natl Cancer Inst (2006) 98(4):262–72. doi: 10.1093/jnci/djj052
24
CharoentongPFinotelloFAngelovaMMayerCEfremovaMRiederDet al. Pan-Cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019
25
KimJYChoiJKJungH. Genome-Wide Methylation Patterns Predict Clinical Benefit of Immunotherapy in Lung Cancer. Clin Epigenet (2020) 12(1):119. doi: 10.1186/s13148-020-00907-4
26
RiazNHavelJJMakarovVDesrichardAUrbaWJSimsJSet al. Tumor and Microenvironment Evolution During Immunotherapy With Nivolumab. Cell (2017) 171(4):934–49.e916. doi: 10.1016/j.cell.2017.09.028
27
MariathasanSTurleySJNicklesDCastiglioniAYuenKWangYet al. Tgfβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501
28
ChenDSMellmanI. Elements of Cancer Immunity and the Cancer-Immune Set Point. Nature (2017) 541(7637):321–30. doi: 10.1038/nature21349
29
PanYLiuGZhouFSuBLiY. DNA Methylation Profiles in Cancer Diagnosis and Therapeutics. Clin Exp Med (2018) 18(1):1–14. doi: 10.1007/s10238-017-0467-0
30
EstécioMRIssaJP. Dissecting DNA Hypermethylation in Cancer. FEBS Lett (2011) 585(13):2078–86. doi: 10.1016/j.febslet.2010.12.001
31
GajewskiTFSchreiberHFuYX. Innate and Adaptive Immune Cells in the Tumor Microenvironment. Nat Immunol (2013) 14(10):1014–22. doi: 10.1038/ni.2703
32
GalonJBruniD. Approaches to Treat Immune Hot, Altered and Cold Tumours With Combination Immunotherapies. Nat Rev Drug Discovery (2019) 18(3):197–218. doi: 10.1038/s41573-018-0007-y
33
FridmanWHZitvogelLSautès-FridmanCKroemerG. The Immune Contexture in Cancer Prognosis and Treatment. Nat Rev Clin Oncol (2017) 14(12):717–34. doi: 10.1038/nrclinonc.2017.101
34
HegdePSKaranikasVEversS. The Where, the When, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition. Clin Cancer Res (2016) 22(8):1865–74. doi: 10.1158/1078-0432.Ccr-15-1507
35
GajewskiTF. The Next Hurdle in Cancer Immunotherapy: Overcoming the Non-T-Cell-Inflamed Tumor Microenvironment. Semin Oncol (2015) 42(4):663–71. doi: 10.1053/j.seminoncol.2015.05.011
36
BatlleEMassaguéJ. Transforming Growth Factor-β Signaling in Immunity and Cancer. Immunity (2019) 50(4):924–40. doi: 10.1016/j.immuni.2019.03.024
37
PeranzoniELemoineJVimeuxLFeuilletVBarrinSKantari-MimounCet al. Macrophages Impede CD8 T Cells From Reaching Tumor Cells and Limit the Efficacy of Anti-PD-1 Treatment. Proc Natl Acad Sci USA (2018) 115(17):E4041–50. doi: 10.1073/pnas.1720948115
38
SunRLimkinEJVakalopoulouMDercleLChampiatSHanSRet al. A Radiomics Approach to Assess Tumour-Infiltrating CD8 Cells and Response to Anti-PD-1 or Anti-PD-L1 Immunotherapy: An Imaging Biomarker, Retrospective Multicohort Study. Lancet Oncol (2018) 19(9):1180–91. doi: 10.1016/s1470-2045(18)30413-3
39
XiaoYMaDZhaoSSuoCShiJXueMZet al. Multi-Omics Profiling Reveals Distinct Microenvironment Characterization and Suggests Immune Escape Mechanisms of Triple-Negative Breast Cancer. Clin Cancer Res (2019) 25(16):5002–14. doi: 10.1158/1078-0432.Ccr-18-3524
40
de CubasAADunkerWZaninovichAHongoRABhatiaAPandaAet al. DNA Hypomethylation Promotes Transposable Element Expression and Activation of Immune Signaling in Renal Cell Cancer. JCI Insight (2020) 5(11):e137569. doi: 10.1172/jci.insight.137569
41
WuXNSuDMeiYDXuMQZhangHWangZYet al. Identified Lung Adenocarcinoma Metabolic Phenotypes and Their Association With Tumor Immune Microenvironment. Cancer Immunol Immunother (2021). doi: 10.1007/s00262-021-02896-6
42
XuJZGongCXieZFZhaoH. Development of an Oncogenic Driver Alteration Associated Immune-Related Prognostic Model for Stage I-II Lung Adenocarcinoma. Front Oncol (2020) 10:593022. doi: 10.3389/fonc.2020.593022
43
ZhaoJXiaoZLiTChenHYuanYWangYAet al. Stromal Modulation Reverses Primary Resistance to Immune Checkpoint Blockade in Pancreatic Cancer. ACS Nano (2018) 12(10):9881–93. doi: 10.1021/acsnano.8b02481
44
FeigCJonesJOKramanMWellsRJDeonarineAChanDSet al. Targeting CXCL12 From FAP-Expressing Carcinoma-Associated Fibroblasts Synergizes With Anti-PD-L1 Immunotherapy in Pancreatic Cancer. Proc Natl Acad Sci USA (2013) 110(50):20212–7. doi: 10.1073/pnas.1320318110
45
GaneshKMassaguéJ. TGF-β Inhibition and Immunotherapy: Checkmate. Immunity (2018) 48(4):626–8. doi: 10.1016/j.immuni.2018.03.037
46
TaurielloDVFPalomo-PonceSStorkDBerenguer-LlergoABadia-RamentolJIglesiasMet al. Tgfβ Drives Immune Evasion in Genetically Reconstituted Colon Cancer Metastasis. Nature (2018) 554(7693):538–43. doi: 10.1038/nature25492
47
LeeSYJeongEKJuMKJeonHMKimMYKimCHet al. Induction of Metastasis, Cancer Stem Cell Phenotype, and Oncogenic Metabolism in Cancer Cells by Ionizing Radiation. Mol Cancer (2017) 16(1):10. doi: 10.1186/s12943-016-0577-4
48
BelleiBMiglianoE. Picardo M. A Framework of Major Tumor-Promoting Signal Transduction Pathways Implicated in Melanoma-Fibroblast Dialogue. Cancers (Basel) (2020) 12(11):3400. doi: 10.3390/cancers12113400
49
LukeJJFlahertyKTRibasALongGV. Targeted Agents and Immunotherapies: Optimizing Outcomes in Melanoma. Nat Rev Clin Oncol (2017) 14(8):463–82. doi: 10.1038/nrclinonc.2017.43
50
SkoulidisFGoldbergMEGreenawaltDMHellmannMDAwadMMGainorJFet al. STK11/LKB1 Mutations and PD-1 Inhibitor Resistance in KRAS-Mutant Lung Adenocarcinoma. Cancer Discovery (2018) 8(7):822–35. doi: 10.1158/2159-8290.Cd-18-0099
51
DongZYZhongWZZhangXCSuJXieZLiuSYet al. Potential Predictive Value of TP53 and KRAS Mutation Status for Response to PD-1 Blockade Immunotherapy in Lung Adenocarcinoma. Clin Cancer Res (2017) 23(12):3012–24. doi: 10.1158/1078-0432.Ccr-16-2554
52
ChengXYinHFuJChenCAnJGuanJet al. Aggregate Analysis Based on TCGA: TTN Missense Mutation Correlates With Favorable Prognosis in Lung Squamous Cell Carcinoma. J Cancer Res Clin Oncol (2019) 145(4):1027–35. doi: 10.1007/s00432-019-02861-y
53
ZhangBWuQLiBWangDWangLZhouYL. M(6)A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol Cancer (2020) 19(1):53. doi: 10.1186/s12943-020-01170-0
54
DiNardoARRajapaksheKNishiguchiTGrimmSLMtetwaGDlaminiQet al. DNA Hypermethylation During Tuberculosis Dampens Host Immune Responsiveness. J Clin Invest (2020) 130(6):3113–23. doi: 10.1172/jci134622
Summary
Keywords
lung adenocarcinoma (LUAD), DNA methylation, tumor microenvironment, immunotherapy, mutation burden
Citation
Yuan D, Wei Z, Wang Y, Cheng F, Zeng Y, Yang L, Zhang S, Li J and Tang R (2021) DNA Methylation Regulator-Meditated Modification Patterns Define the Distinct Tumor Microenvironment in Lung Adenocarcinoma. Front. Oncol. 11:734873. doi: 10.3389/fonc.2021.734873
Received
01 July 2021
Accepted
16 August 2021
Published
06 September 2021
Volume
11 - 2021
Edited by
Giannis Mountzios, National and Kapodistrian University of Athens, Greece
Reviewed by
Petros Christopoulos, Heidelberg University Hospital, Germany; Jules Louis Derks, Maastricht University Medical Centre, Netherlands
Updates
Copyright
© 2021 Yuan, Wei, Wang, Cheng, Zeng, Yang, Zhang, Li and Tang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Renkuan Tang, 100456@cqmu.edu.cn
†These authors have contributed equally to this work
This article was submitted to Thoracic Oncology, a section of the journal Frontiers in Oncology
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.