Molecular Characterization of the Clinical and Tumor Immune Microenvironment Signature of 5-methylcytosine-Related Regulators in non-small Cell Lung Cancer
Taisheng Liu1† Liyi Guo2† Guihong Liu3† Xiaoshan Hu4† Xiaoning Li1 Jinye Zhang1 
Zili Dai5 Peng Yu5 Ming Jiang6 Jian Wang1* 
Jian Zhang5*- 1Department of Thoracic Surgery, Affiliated Cancer Hospital and Institute of Guangzhou Medical University, Guangzhou, China
- 2The Sixth People’s Hospital of Huizhou City, Huiyang Hospital Affiliated to Southern Medical University, Huizhou, China
- 3Department of Radiation Oncology, DongGuan Tungwah Hospital, Dongguan, China
- 4Department of Internal Medicine of Oncology, Affiliated Cancer Hospital and Institute of Guangzhou Medical University, Guangzhou, China
- 5Department of Radiation Oncology, Affiliated Cancer Hospital and Institute of Guangzhou Medical University, State Key Laboratory of Respiratory Diseases, Guangzhou Institute of Respiratory Disease, Guangzhou, China
- 6Department of Breast Surgery, Affiliated Cancer Hospital and Institute of Guangzhou Medical University, Guangzhou, China
Background: DNA methylation is an important epigenetic modification, among which 5-methylcytosine methylation (5mC) is generally associated with tumorigenesis. Nonetheless, the potential roles of 5mC regulators in the tumor microenvironment (TME) remain unclear.
Methods: The 5mC modification patterns of 1,374 lung adenocarcinoma samples were analyzed systematically. The correlation between the 5mC modification and tumor microenvironment cell infiltration was further assessed. The 5mCscore was developed to evaluate tumor mutation burden, immune check-point inhibitor response, and the clinical prognosis of individual tumors.
Results: Three 5mC modification patterns were established based on the clinical characteristics of 21 5mC regulators. According to the differential expression of 5mC regulators, three distinct 5mC gene cluster were also identified, which showed distinct TME immune cell infiltration patterns and clinical prognoses. The 5mCscore was constructed to evaluate the tumor mutation burden, immune check-point inhibitor response, and prognosis characteristics. We found that patients with a low 5mCscore had significant immune cell infiltration and increased clinical benefit.
Conclusion: This study indicated that the 5mC modification is involved in regulating TME infiltration remodeling. Targeting 5mC modification regulators might be a novel strategy to treat lung cancer.
Introduction
Lung cancer is the primary cause of cancer-related deaths worldwide (Siegel et al., 2020) (NSCLC), accounting approximately for 85% of newly diagnosed lung cancer cases, is classified into lung adenocarcinoma (LUAD) and lung squamous carcinoma (LUSC) (Curran et al., 2011). For unresectable advanced NSCLC, a combination of radiotherapy and chemotherapy has been the most common first-line treatment (Yoda et al., 2019), and impressive clinical success has been observed using targeted therapies (Treat, 2005; Yuan et al., 2019; Alexander et al., 2020). Unfortunately, most NSCLC patients will suffer the relapse within 1 year (Fountzilas et al., 2021). Thus, understanding the mechanism and identifying novel targets to treat NSCLC remain an urgent clinical need.
Immunotherapies represent a promising advance in cancer treatment (Lussier et al., 2021). The immune checkpoint inhibitors (ICI), including programmed death-ligand 1 (PD-L1), programmed cell death 1 (PD-1), and cytotoxic T-lymphocyte antigen-4 (CTLA-4), combined with chemoradiotherapy, have been approved or are being widely evaluated in clinical trials (Grant et al., 2021). However, targeting PD-1 or PD-L1 has demonstrated durable efficacy only in a subset of patients with NSCLC (Jazieh et al., 2021). Thus, it is important to determine the underlying mechanisms with the aim of improving the curative effect.
DNA methylation is an epigenetic modification that is associated with regulating cell differentiation and tissue development (Smith and Meissner, 2013; Slieker et al., 2015). Dysregulation of DNA methylation patterns are important characteristics of several diseases, including cancers (Li et al., 2013; Božić et al., 2021; Cristall et al., 2021; Miyakuni et al., 2021). 5-Methylcytosine (5mC), a type of DNA methylation, was firstly reported by Wyatt, (1951). DNA 5mC methylation is the classic epigenetic process, which is controlled by “writers” (DNA methyltransferases), “erasers” (DNA methyltransferases), and “readers” (Ito et al., 2011; Du et al., 2015; Lio et al., 2020). With the discovery of 5mC regulators, recent studies suggested that DNA cytosine modifications may act as epigenetic markers in tumorigenesis (Wu and Zhang, 2010; Cavalcante et al., 2020; Jiang, 2020; Mo et al., 2020) and can regulate tumor microenvironment (TME) infiltrating cells (Chen et al., 2020; Zhao et al., 2021; Onodera et al., 2021). However, the comprehensive roles of TME cell infiltration directed by 5mC regulators in NSCLC remain unclear.
In this study, we evaluated 5mC methylation patterns comprehensively by analyzing genomic information of 1374 LUAD samples, and correlated the 5mC methylation pattern with the characteristics of TME cell infiltration. We identified three 5mC methylation patterns, and revealed that 5mC methylation mediation of TME cell infiltration characteristics was closely associated with the immune response phenotype, indicating the 5mC methylation played an important role in modifying TME characteristics. Furthermore, the 5mCscore could be applied as a promising biomarker to predict immune response and clinical outcome in NSCLC.
Materials and Methods
Dataset Acquisition and Processing
Supplementary Figure S1 shows the workflow of the this study. mRNA expression with clinical and survival information were downloaded from Gene Expression Omnibus (GEO) and GDC data portal. Patients without clinical survival information were excluded. Five eligible lung adenocarcinoma cohorts (GSE19188, GSE31210, GSE37745, GSE50081, and TCGA-LUAD [lung adenocarcinoma data from The Cancer genome Atlas (TGCA)]) were included for further analysis (Supplementary Table S1). For background correction and normalization, the Robust Multichip Average algorithm was used to uniformly process the raw. CEL files of the four GEO datasets (Gautier et al., 2004). Next, a GEO meta-cohort were created by merging the GEO datasets using the R sva package (Leek et al., 2012).
Twenty-one 5mc regulators, including three writers (DNA methyltransferase 1 (DNMT1), DNA methyltransferase 3 Alpha (DNMT3A), DNA methyltransferase 3 beta (DNMT3B)), three erasers (tet methylcytosine dioxygenase 1 (TET1), tet methylcytosine dioxygenase 2 (TET2), tet methylcytosine dioxygenase 3 (TET3)), and 15 readers (methyl-cpg binding domain protein 1 (MBD1), methyl-cpg binding domain protein 2 (MBD2), methyl-cpg binding domain protein 3 (MBD3), methyl-cpg binding domain protein 4 (MBD4), methyl-cpg binding protein 2 (MECP2), nei like dna glycosylase 1 (NEIL1), nth like dna glycosylase 1 (NTHL1), single-strand-selective monofunctional uracil-dna glycosylase 1 (SMUG1), thymine dna glycosylase (TDG), ubiquitin like with phd and ring finger domains 1 (UHRF1), ubiquitin like with phd and ring finger domains 2 (UHRF2), uracil dna glycosylase (UNG), zinc finger and btb domain containing 33 (ZBTB33), zinc finger and btb domain containing 34 (ZBTB34), zinc finger and btb domain containing 4 (ZBTB4)) (Chen et al., 2020), and 23 tumor immune related cells from published studies (Zhang et al., 2020a; Zhao et al., 2021), were included for analysis. The transcriptomics data, single nucleotide variant (SNV), copy number variation (CNV), and 5mC phenotypic data were collected using the UCSC Xena database (https://xenabrowser.net) and the GDC data portal.
Unsupervised Clustering of 21 5mC Regulators
To identify 5mC regulator-mediated modification sub-clusters, unsupervised consensus clustering was used to cluster tumor samples into sub-clusters based on the expression levels of the 21 5mC regulators. To ensure the stability of the clusters, the parameters of clustering were as follows: number of repetitions = 1,000 bootstraps, clustering algorithm = k-means method, pFeature = 1.0, pItem = 0.8. The cluster with the most significant survival difference was included for further analysis.
Gene Set Variation Analysis and Functional Annotation
To explore the biological behavior among the different 5mC modification patterns, their pathway scores were evaluated using gene set variation analysis (GSVA) using the R GSVA package (Hänzelmann et al., 2013), with the “c2. cp.kegg.v7.4. symbols” gene set as the background. Differential pathways were further screened using p < 0.05 in the R package limma.
Estimation of the Tumor Microenvironment
To identify the TME cell infiltration in LUAD samples, the relative abundances of immune cells were quantified using the single-sample gene-set enrichment analysis (ssGSEA) algorithm. According to the method revealed by Charoentong et al. (2017b), various kinds of immune cells, including regulatory T cells, activated CD8+ T cells, dendritic cells, and B cells, were evaluated. The relative abundance of TME infiltrating cells in clinical samples was represented by the enrichment scores.
Differentially Expressed Genes
To identify 5mC-related differentially expressed genes (DEGs), based on the expression levels of 21 5mC regulators, three distinct 5mC modification patterns were identified in the patients with LUAD. The empirical Bayesian approach of R package limma package was used for the difference analysis (Ritchie et al., 2015), which screened out 324 DEGs, 246 DEGs and 144 DEGs according to p < 0.001, p < 0.0005 and p < 0.0001. p < 0.0005 was most suitable for subsequent analysis.
Construction of 5mC Gene Signatures
Considering the heterogeneity and complexity of tumors, and according to the method used by Zhang J. et al. (2020), the 5mCscore was developed to quantify the modification pattern of individual patients with LUAD based on the identified DEGs. A univariate Cox regression model was used for the prognostic analysis of each gene in the 5mC signatures. We obtained 103 genes related to prognosis from among the 246 DEGs, and then principal component analysis (PCA) was performed, scored as PCi1 and PCi2. This approach had advantage of focusing the score on the set with the largest block of well correlated (or anticorrelated) genes in the set, while down-weighting contributions from genes that do not track with other set members. The 5mC score of each patient was calculated as follows:
Evaluation of Immune-Checkpoint Inhibitor Genomic and Clinical Information
To explore the application of the 5mC score to predict immune-checkpoint inhibitor (ICI) efficacy, the expression data and clinical annotations of the immunotherapeutic cohort of atezolizumab (IMvigor210 cohort) were downloaded from the website based on the Creative Commons 3.0 License (http://research-pub.Gene.com/imvigor210corebiologies) (Mariathasan et al., 2018).
Statistical Analysis
Correlation coefficients between the expression of 5mC regulators and the TME immune infiltration cells was conducted using the Spearman method and distance correlation analysis. The Wilcoxon test was used to analyze the difference between two groups. The Kruskal–Wallis test and one-way analysis of variance (ANOVA) were used to analyze difference among three or more groups. The log-rank test and the Kaplan–Meier (KM) method were applied to evaluate the survival time. A statistical two-sided p value < 0.05 was considered as having significance. All data processing in this study was done using R 3.6.1 software.
Results
Genetic Variation and Expression Analysis of 5mC Methylation Regulators
According to the map described in Figure 1A, in this study, 21 5 mC methylation regulators (writers: DNMT1, DNMT3A and DNMT3B; erasers: TET1, TET2, TET3; readers: MBD1, MBD2, MBD3, MBD4, MECP2, NEIL1, NTHL1, SMUG1, TDG, UHRF1, UHRF2, UNG, ZBTB33, ZBTB38, and ZBTB4) were identified (Supplementary Table S2). To determine genetic alternations, we firstly evaluated the SNV variation frequency of the genes encoding the 21 5mC methylation regulators. As shown in Figure 1B, Among 561 LUAD samples, 21.39% of 5mC regulators had mutations. The main mutation type was missense_mutation. However, the mutation frequency of individual regulators only ranged from 0 to 4%. The CNV frequency of the 5 mC regulators showed that MECP2, SMUG1, DNMT3B, ZBTB33, and NTHL1 had distinct CNV amplification, with frequencies of 11.71, 6.13, 5.58, 4.68, and 4.68%, respectively. MBD3, UHRF1, MBD1, UHRF2, and TDG had a CNV deletion, with frequencies of 6.30, 5.40, 5.58, 5.23, and 4.86%, respectively (Figure 1C and Supplementary Table S3). The distribution analysis of CNV alterations on 23 chromosomes showed that their distribution among the 21 5mC regulators was scattered and unorganized (Figure 1D). Survival analysis indicated that high expression of DNMT3B, MDB2, MDB3, SMUG1, TDG, HURF1, UNG, and ZBTB38 were associated with poor survival of LUAD (p < 0.05); while, high expression of MDB4, MECP2, NEIL1, TET2, UHRF2, and ZBTB4 were associated with better survival of LUAD (p < 0.05, Supplementary Figure S2 and Supplementary Table S4).
FIGURE 1. Genetic landscape and expression analysis of 5mC regulators in LUAD. (A) Schematic diagram of 5mC DNA methylation mediated by 21 5mC regulators. (B) The mutation frequency of 21 5mC regulators in the TCGA-LUAD cohort. The column indicates individual patients. The upper barplot shows the TMB, The number on the right indicates the mutation frequency. The right barplot shows the percentage of mutation type in each regulator. The stacked barplot shows the fraction of conversions. (C) The CNV variation frequency of 5mC regulators in the TCGA-LUAD cohort. The height of the column indicates the alteration frequency of the regulators. The green dot is the deletion frequency; The red dot is the amplification frequency. (D) The locations of CNV alterations of 5mC regulators in the TCGA-LUAD cohort. CNV, copy number variation; 5mC, 5-methylcytosine; LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; TMB, tumor mutation burden.
Identification of 5mC Methylation-Related Phenotypes
To determine the roles of interaction among 5 mC methylation regulators in LUAD, correlation analysis among the 21 5 mC regulators was performed, which showed that there was a strong positive correlation among most of the regulators (Supplementary Figure S3A and Supplementary Table S5). The prognostic values of the 21 5 mC regulators in LUAD were evaluated using a univariate Cox regression model (Supplementary Figure S3B). As shown in Figure 2A, MDB4, MECP2, NEIL1, TET2, ZBTB4, and ZBTB33 were favorable factors for overall survival (OS), while DNMT1, DNMT3A, DNMT3B, TET1, TET3, SMUG1, TDG, UHRF1, UHRF2, UNG, and ZBTB38 were risk factors for OS. Significant negative correlations were obtained for UHRF1 and DNMT1, TDG and DNMT3A, TDG and UNG, MECP2 and ZBTB33, MECP2 and TET2, and TET2 and UHRF2 (p < 0.001). On the other hand, several erasers and readers also showed significant negative correlations: NTHL1 and TET2, NTHL1 and TET3, and MBD3 and TET2 (p < 0.001) (Supplementary Tables S6–7). Using unsupervised clustering analysis, three distinct 5mC modification patterns were identified based on the expression of 21 5mC regulators (Supplementary Figure S4). Prognostic analysis of the three 5mC modification clusters revealed a particularly prominent survival advantage for the 5mC cluster-B modification pattern (Figure 2B and Supplementary Table S8; p = 0.001). The results showed that cross-talk among the 5mC modification regulators might be involved in the formation of the 5mC modification and in the characteristics of TME cell infiltration.
FIGURE 2. Prognostic and biological characteristics of 5mC modification patterns in LUAD. (A) The interaction among 5mC regulators in the TCGA-LUAD cohort. The circle size indicates the effect of each regulator on prognosis. Green dots indicate favorable factors for prognosis; Purple dots indicate risk factor for prognosis. The lines linking regulators indicate their interactions, and thickness show the correlation strength between the regulators. Negative correlations are marked in blue and positive correlation in red. (B) Survival analyses of 5mC modification patterns in the TCGA-LUAD cohort, including 500 cases in 5mC cluster A (n = 139), 5mC cluster B (n = 187), and 5mC cluster C (n = 174) (p < 0.0001, Log-rank test). (C) Cluster analysis of 21 5mC regulators among the three 5mC modification patterns. (D) Gene ontology (GO) analysis of 21 5mC regulators among the three 5mC modification patterns.
Tumor Microenvironment Cell Infiltration Characteristics in the 5mC Methylation Clusters
To identify the potential function of the differentially expressed 5mC regulators, cluster analysis was first performed. As shown in Figure 2C, the 21 5mC regulators had a distinct distribution among the three 5mC clusters. Gene ontology (GO) analysis was performed to identify the biological process (BP), cellular component (CC), and molecular function (MF) of the regulators. The aberrantly expressed 5mC regulators were mainly enriched for GO terms related to regulation of mitotic nuclear division, chromosome segregation, and nuclear division (BP); chromosomal region, condensed chromosome/centromeric region, and kinetochore (CC); and ATPase activity, DNA helicase activity, and helicase activity (MF) (Figure 2D). To further identify the potential behaviors, GSVA enrichment analysis was performed, as shown in Supplementary Table S9. To further exlpore unsupervised consensus clustering of all tumor samples for the molecular classification of LUAD. The optimal number of clusters was determined by the K value. After assessing relative changes in the area under the cumulative distribution function curve and consensus matrix heatmap, we selected a three-cluster solution (K = 3), which showed no appreciable increase in the area under the cumulative distribution function curve (Supplementary Figure S4). 5mC cluster A was markedly enriched in damage repair-related pathways, such as base excision repair, DNA replication, spliceosome, and RNA polymerase. 5mC cluster B was prominently related to immune activation-related pathways, such as the JAK-STAT signaling pathway, the T cell receptor signaling pathway, and the calcium signaling pathway. 5mC cluster C was mostly associated with carcinogenic activation and damage repair pathways, such as, the p53 signaling pathway, basal transcription factors, spliceosome, RNA degradation, DNA replication, base excision repair, homologous recombination, DNA replication, and mismatch repair (Figures 3A–C). Based on the expression levels of these 21 5mC regulators, the three 5mC modification patterns could be partially differentiated using PCA (Figure 3D). TME cell infiltration analysis showed 5mC cluster B was associated with activated B cells, activated dendritic cells, mast cells, natural killer T cells, and neutrophils (Figure 3E and Supplementary Table S10, p < 0.001). 5mC cluster C was remarkably rich in immune cell infiltration including myeloid-derived suppressor cells (MDSCs), regulatory T cells, type 1 T helper cells, type 2 T helper cells, and type 17 T helper cell (Figure 3E and Supplementary Table S10, p < 0.001). Prognosis analysis showed that patients with different 5mC modification patterns also had a matching survival advantage (Figure 2B, p = 0.001). Based on the above results, cluster A, characterized by innate immune cell infiltration, was defined as an immune-excluded phenotype; cluster B, characterized by adaptive immune cell infiltration and immune activation, was defined as an immune-inflamed phenotype; and cluster C, characterized by the inhibition of immunity, was defined as an immune-desert phenotype.
FIGURE 3. GSVA enrichment analysis and TME cell infiltration characteristics of 5mC modification patterns. (A–C) The states of biological pathways among the three 5mC modification patterns enriched by GSVA analysis The general biological processes are shown as a heatmap, red represents activated pathways and blue represents inactivated pathways. (A) 5 mC cluster A vs 5mC cluster B; (B) 5mC cluster B vs 5mC cluster C; (B) 5 mC cluster A vs 5mC cluster C. (D) Principal component analysis of the 5mC modification patterns. (E) The abundance of TME infiltrating cells in the 5mC modification patterns (*p < 0.05; **p < 0.01; ***p < 0.001; ns means not significant).
Identification of 5mC Methylation Gene Signature
To further identify the potential function of each m5C modification pattern, we determined 246 m5C phenotype-related DEGs (Supplementary Table S11). GO analysis showed that the 246 DEGs were associated with cell cycle, RNA transport, spliceosome, DNA replication, base excision, and human T-cell leukemia virus 1 infection (Supplementary Figure S5A and Supplementary Table S12). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated that the 5mC gene clusters were involved in DNA transcription and translation (Supplementary Figure S5B and Supplementary Table S13). To further determine the potential regulation mechanism, unsupervised clustering analyses was performed to identify the genomic subtypes based on the 103 prognostic genes from the 246 5mC phenotype-related DEGs. The results showed that three distinct 5 mC genomic phenotypes (5mC gene Cluster A–C) could be identified (Figure 4A and Supplementary Table S5C–J). These results indicated that the 5mC methylation modification patterns did exist in LUAD and three distinct 5mC gene clusters were characterized by different signature genes. Cluster analysis showed that 178 of 504 patients with LUAD were clustered in 5mC gene cluster C, which was associated with better prognosis. Patients with LUAD with 5mC gene cluster B (n = 135) had poorer prognosis. 5mC gene cluster A, with 191 patients clustered, had an intermediate prognosis (Figure 4B, p < 0.001). The expression levels of the 5mC regulators among the 5mC gene clusters were distinctly different (Figure 4C).
FIGURE 4. Construction of 5mC gene signatures. (A) Unsupervised clustering of overlapping 5mC phenotype-related genes in the TCGA-LUAD cohort to classify patients into different genomic subtypes, termed as 5mC gene cluster (A–C), respectively. The gene clusters, 5mC clusters, tumor stage, survival status, sex, and age were used as patient annotations. (B) Overall survival of patients with the three 5mC modification genomic clusters in the TCGA-LUAD cohort, including 504 cases in 5mC gene cluster A (n = 191), 5mC gene cluster B (n = 135), and 5mC gene cluster C (n = 17) (p < 0.0001, Log-rank test). (C) The expression of 21 5mC regulators in the three gene clusters (*p < 0.05; **p < 0.01; ***p < 0.001; ns means not significant). (D) Alluvial diagram showing the changes in 5mC clusters, 5mC gene cluster, 5mCscore, and survival. (E) Correlations between the 5mCscore and the known gene signatures in the TCGA-LUAD cohort using Spearman analysis. Negative correlations are marked in blue and positive correlation in red. (F) Differences in the 5mCscore among three gene clusters in the TCGA-LUAD cohort (***p < 0.001, Kruskal-Wallis test). (G) Differences in the 5mCscore among three the 5mC modification patterns in the TCGA-LUAD cohort (***p < 0.001, Kruskal-Wallis test).
Clinical Characteristics of 5mCscore Phenotypes
To better explore the pattern of 5mC modification in individual patients, based on the 5mC phenotype-related genes (Supplementary Table S14), the 5mCscore was used to quantify the 5mC modification patterns of individual patients with LUAD. An alluvial diagram was applied to clarify the attributed changes of the LUAD patients. As shown in Figure 4D, the 5mC modification patterns clusters were almost consistent with the 5mC gene clusters, i.e., the 5mC gene cluster B group patients mainly had a low 5mCscore, which was associated with poor survival. To determine the roles of 5mC-related phenotypes in immune regulation, correlation analysis showed that the 5mCscore was associated positively with most TME infiltrating cells (Figure 4E). The Kruskal–Wallis test revealed there was a significant difference in the 5mCscore among the 5mC gene clusters. 5mC gene cluster C showed the highest median 5mCscore, while 5mC gene cluster B had the lowest median 5mCscore, which indicated that a high 5mCscore was closely associated with immune activation-related signatures, whereas a low 5mCscore was associated with immune inactivation-related signatures (Figure 4F, p < 0.001). More importantly, compared with the other clusters, 5mC modification cluster C presented the lowest median 5mCscore, and 5mC modification cluster B showed the highest 5mCscore (Figure 4G, p < 0.001). These results indicated that a high 5mCscore correlated significantly with immune-activation and the 5mCscore could be used to identify the 5mC modification patterns in LUAD, and further assess the characteristics of TME cell infiltration of individual tumors.
To further validate the value of the 5mCscore, patients in the TCGA cohort were divided into low or high 5mCscore groups. Prognosis analysis showed that patients with a high 5mCscore showed a better survival benefit (Figure 5A, p < 0.001). Four GEO datasets (GSE19188, GSE31210, GSE37745, and GSE50081, Supplementary Table S1) were integrated into one meta-cohort. Survival analysis in the GEO meta-cohort also identified that a high 5mCscore was linked to a better clinical outcome (Figure 5B, p < 0.001). These results indicated that the 5mCscore could act as an independent prognostic biomarker to evaluate patient outcomes. To explore the effect of clinical characteristics on the 5mCscore, the subgroups of clinical characteristics were further analyzed. A significant distribution difference of a high 5mCscore was observed for gender (59% in female vs 41% in male, p = 0.0054; Figures 5C,D), smoking status (41% vs 67% for ever smoking, P = 5e-05; Figures 5G,H), stage I–II (83% vs 54% for stage I, p = 3.8e-08; Figures 5I,J), and genetic mutations (63% vs 41% for EGFR mutations, p = 0.00019; 25% vs 42% in EGFR/KRAS/ALK mutations, p < 0.001; Figures 5K-L). However, there were no 5mCscore differences between age (≤65) and age (>65) (Figures 5E,F, p = 0.6). To assess the value of clinical characteristics, patients in the TCGA-LUAD cohort were further stratified by age (≤65/> 65), sex (female/male), T stage (T1–2/T3–4), N stage (N0–1/N2–3), M stage (M0/M1), and clinical stage (I–II/III–IV). We found that the clinical characteristics, particularly T1–2, N0–1, M0, and I–II clinical stages, could be clearly divided into high- and low-risk subgroups (Supplementary Figure S6). These results indicated that multiple clinical characteristics can have an effect on the 5mCscore, which led to the heterogeneity of 5mC regulators in LUAD.
FIGURE 5. Prognostic and genetic characteristics between high and low 5mCscore groups. (A) Survival analysis of the 5mCscore in the TCGA-LUAD cohort (p < 0.0001, Log-rank test). (B) Survival analysis of the 5mCscore in the GEO-meta cohort (p < 0.0001, Log-rank test). (C) Sex proportion between the high- and low-5mCscore groups. (D) The 5mCscore difference between females and males. (E) Age proportion between the high- and low-5mCscore groups. (F) The 5mCscore difference between age (≤65) and age (>65). (G) Smoking status proportion between the high- and low-5mCscore groups. (H) The 5mCscore difference between smoking status (ever) and smoking status (never). (I) Clinical stage status proportion between the high- and low-5mCscore groups. (J) The 5mCscore difference between stage I and stage II. (K) The 5mCscore difference between genetic mutations (−) and genetic mutations (+). (L) Genetic mutations status proportion between high- and low-5mCscore groups.
The Potential of the 5mCscore to Predict the Response to anti-PD-L1 Immunotherapy
The above analyses demonstrated the impact of 5mCscore regulators on the TME, as well as on the prognosis in patients with LUAD. The genetic characteristics of the patients in different 5mCscore groups were further explored. As shown in Figures 6A,B and Supplementary Table S15, the somatic mutation landscapes in the high and low 5mCscore groups had a distinct difference. The mutation frequency was 77.35% in the high 5mCscore group and 94.89% in the low 5mCscore group. Specifically, except for KRAS, TP53 (18% vs 58%), TTN (20% vs 53%), MUC16 (28% vs 45%), and RYR2 (22% vs 40%) had important differences between the high and low 5mCscore groups (Figure 6B). Besides, patients with a low 5mCscore showed a significantly higher tumor mutation burden (TMB) and PD-L1 expression than patients with a high 5mCscore (Figures 6C,D and Supplementary Table S16). 5mC gene cluster C showed lower PD-L1 expression and a lower TMB than 5mC gene cluster B. Correlation analysis further identified that the TMB and PD-L1 expression were related negatively with the 5 mCscore (Figures 6E,F, p < 0.001). These results revealed a significant association between the 5mCscore and the TMB and PD-L1 expression. These factors are important parameters in the assessment of immunotherapy outcomes. However, the survival analysis associated with the TMB found that there was no difference between the high and low TMB groups (Figure 6G, p = 0.082). Next, the crosstalk between the 5mCscore and TMB in terms of patient survival was investigated. The high 5mCscore and high TMB group had better survival than the low 5mCscore and high TMB group. The low 5mCscore and low TMB group was associated with poorer survival relative to those with a high 5mCscore and low TMB (Figure 6H, p < 0.001).
FIGURE 6. Characteristics of 5mCscore in the TCGA molecular subtypes and tumor somatic mutations. (A,B) Waterfall plot of tumor somatic mutations established by those with a high 5mCscore (A) and a low 5mCscore (B). Each column represents individual patients. The upper barplot shows the TMB, the number on the right indicates the mutation frequency in each gene. The right barplot shows the proportion of each variant type. (C) Tumor somatic mutation between high 5mCscore and low 5mCscore groups. (D) PD-L1 expression difference between high 5mCscore and low 5mCscore groups. (E) The correlation analysis between tumor somatic mutation and the 5mCscore. (F) The correlation analysis between PD-L1 expression and the 5mCscore. (G) Survival analysis of tumor somatic mutations in the TCGA-LUAD cohort (p < 0.0001, Log-rank test). (H) Survival analyses for patients stratified by both the 5mCscore and the tumor somatic mutation burden using Kaplan–Meier curves (p < 0.0001, Log-rank test).
To explore the potential roles of the 5mCscore in clinical immune therapy of lung cancer, we investigated whether the 5mCscore could predict patients’ response to PD-L1 (atezolizumab) therapy based on the PD-L1 immunotherapy cohort (IMvigor210). Compared with those with a high 5mCscore, patients with a low 5mCscore had significant therapeutic advantages and clinical responses to anti-PD-L1 immunotherapy (Figures 7A,B and Supplementary Figure S7A−B, p = 0.0015). The low 5mCscore group had a higher immune cells 2 (IC2) score (38% vs 16%) and a lower tumor cells 2+ (TC2+) score (77% vs 96%) than the high 5mCscore group, 5mCscore was significantly associated with the enrollment ICs and suppression of TCs (Figures 7C,D and Supplementary Figure S7C−D). These results identified that the 5mCscore played a non-negligible role in regulating TME immune cell infiltration. We further investigated different immune phenotypes among the high and low 5mCscore groups and found that a higher 5mCscore was markedly associated with exclusion and desert immune phenotypes, in which an antitumor effect is difficult to exert using ICI therapy (Figures 7E,F). Patients with a low 5mCscore exhibited significant clinical benefits and a markedly prolonged survival (Figure 7G, p = 0.015). These results clarified that 5mC modification patterns are significantly associated with immune phenotypes and PD-L1 expression, and that the 5mCscore could be a prominent biomarker to predict the response to ICI therapy.
FIGURE 7. The role of the 5mCscore in anti-PD-L1 immunotherapy. (A) The proportion of patients with a response to ICI in the low or high 5mCscore groups. Responder/Nonresponder: 26%/74% in the low 5mCscore groups and 8%/92% in the high 5mCscore groups. (B) 5mCscore differences between responders and nonresponders. (C) IC infiltration proportion between high 5mCscore and low 5mCscore groups. (D) 5mCscore differences between different IC subgroups. (E) Immune phenotype proportion between high 5mCscore and low 5mCscore groups. (F) 5mCscore differences among the immune-desert phenotype, immune-excluded phenotype, and immune-inflamed phenotype. (G) Survival analyses for low (n = 291) and high (n = 57) 5mCscore patient groups in the anti-PD-L1 immunotherapy cohort using Kaplan–Meier curves (IMvigor210 cohort; p = 0.015, Log-rank test). CR, complete response; IC, immune cell; ICI, immune check-point inhibitor; PD, progressive disease; PR, partial response; SD, stable disease.
Discussion
DNA 5mC methylation is a dynamic and reversible post-transcriptional modification regulated by 5mC related regulators (Mayer et al., 2000; Oswald et al., 2000; Wu et al., 2020). Recent research highlighted the biological importance of 5mC modification on immune cell infiltration and tumor suppression (Schübeler, 2015; Dor and Cedar, 2018; Weng et al., 2021). However, most studies focused only on a single TME cell type or one 5mC related regulator, and the comprehensive roles of 5mC regulators on TME infiltration characteristics are not fully elaborated. Thus, further clarification of the potential roles of 5mC modification patterns in the infiltration of TME cells will raise our awareness of the effects of the heterogeneity and complexity of the TME on the response to ICI therapy and provide a novel biomarkers to evaluate the ICI response and predict prognosis.
Herein, three distinct 5mC methylation modification patterns were identified based on 21 5mC regulators. The patterns had significantly distinct TME cell infiltration characteristics. Based on the identified 246 5mC phenotype-related DEGs, three genomic clusters of 5mC-related genes were further identified, which were also validated for their association with transcription modification and immune infiltration. Recent studies had shown that DNA methylation can be involved in the maintenance and reinforcement of T cell exhaustion gene signatures (Pauken et al., 2016; Gate et al., 2018). In murine antigen-specific CD8 T cells, DNMT3A-mediated methylation impaired T cell expansion and led to immune cell exhaustion under treatment with anti-PD-1 via repression the expression of key genes (Ghoneim et al., 2017). By contrast, in the context of T cell exhaustion, the involvement of DNA methylation in the reprogramming of the T cells has also been reported (Araki et al., 2013), such as demethylation of the PD1 promoter resulting in permanent CD8+T cell exhaustion. Uhrf1-mediated tnf-α gene methylation controlled proinflammatory macrophages in experimental colitis resembling inflammatory bowel disease (Qi et al., 2019). In breast cancer, ZBTB33 subcellular partitioning functionally linked LC3A/B, the tumor microenvironment, and cancer survival (Singhal et al., 2021). These results indicated that the 5mC modification is intimately involved in shaping TME landscapes.
Epigenetic alterations are associated extensively with the immune response and tumore evasion. A DNA methylation signature (the EPIMMUNE signature) has been identified as an epigenetic biomarker of the response to ICI. The multicenter and retrospective analysis revealed that the EPIMMUNE signature could predict the response to anti-PD-1 treatment in non-small-cell lung cancer (Seremet et al., 2016). In metastatic melanoma treated with CTLA-4 blockers, responders and non-responders to ICI had a differential DNA methylation pattern (Chida et al., 2021). To better understand the individual heterogeneity of TME-meditated 5mC modification patterns, the 5mCscore was established to assess the 5mC modification pattern of individuals with LUAD. 5mC gene cluster C, characterized by an immune inflamed phenotype, exhibited a higher 5mCscore, and 5mC gene cluster B, characterized by an immune excluded phenotype, had a lower 5mCscore. These results revealed the 5mCscore was a useful biomarker to comprehensively assess individual tumor 5mC modification patterns, which could be used to evaluate TME immune cell infiltration patterns. Prognosis analyses also identified that the 5mCscore was an independent prognostic biomarker in LUAD.
Alterations in 5mC regulatory genes might also be associated with variations in LUAD. In this study, we identified twenty driver genes, including TP53, TTN, MUC16, RYR2, and CSMD3. Moreover, variations in KRAS were associated significantly with alterations in 5mC regulatory genes. As an oncogene, KRAS mutations were reported frequently in a variety of tumors, including colorectal cancer (Prior et al., 2012), pancreatic cancer (Arner et al., 2019), and bladder cancer (Santha et al., 2020). Recent studies identified that KRAS might have a critical role in the immunoregulation of NSCLC (Li et al., 2021; Wang et al., 2021). Our data also revealed that the 5mCscore had a markedly negative correlation with PD-L1 expression and the TMB. The 5mCscore integrating the TMB could be the more effective biomarker to predict ICI response. We also identified the predictive value of the 5mCscore in the IMvigor210 cohort. The 5mCscore between non-responders and responders was significantly different. These results provided new insights to clarify different tumor immune phenotypes and improve the clinical response to ICI therapy.
Conclusion
In summary, we comprehensively analyzed the potential mechanisms of 5mC methylation modification during the regulation of the TME. 5mC modification patterns contributed to the heterogeneity and complexity of the TME in LUAD, which was significantly associated with TMB, PD-L1 expression, and immune phenotypes. 5mCscore could act as a biomarker to predict a patient’s response to ICI therapy.
Data Availability Statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author Contributions
TL and JZ performed all experiments, prepared figures, and drafted the manuscript. LG, XH, XL, and JZ participated in data analysis and interpretation of the results. GL, JZ, XH, ZD, PY, MJ, and JW designed the study and participated in the data analysis. All authors have read and approved the manuscript.
Funding
This study was supported by grants from the National Natural Science Foundation of China (82003212), Discipline Construction Project of Guangzhou Medical University during the 14th Five-Year Plan (06-410-2107181), Guangzhou Key Medical Discipline Construction Project Fund (02-412-B205002-1004042), and the Medical and Health Technology Projects of Guangzhou (2015A011086).
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.
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.
Acknowledgments
We thank the staff members of The Cancer Genome Atlas for their involvement in the cBioPortal for Cancer Genomics Program.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.779367/full#supplementary-material
Abbreviations
5mC: 5-methylcytosine methylation; TME: tumor microenvironment; NSCLC: (Siegel et al., 2020); LUAD: lung adenocarcinoma; LUSC: lung squamous carcinoma; ICI: The immune checkpoint inhibitors; PD-L1: programmed death-ligand 1; PD-1: programmed cell death 1; CTLA-4: cytotoxic T-lymphocyte antigen-4; GEO: Gene Expression Omnibus; SNV: single nucleotide variant; CNV: copy number variation; GSVA: gene set variation analysis; ssGSEA: single-sample gene-set enrichment analysis; DEGs: differentially expressed genes; ANOVA: one-way analysis of variance; KM: Kaplan–Meier; TMB: tumor mutation burden; GO: Gene ontology; BP: biological process; CC: cellular component; MF: molecular function; MDSCs: myeloid-derived suppressor cells.
References
Alexander, M., Kim, S. Y., and Cheng, H. (2020). Update 2020: Management of Non-small Cell Lung Cancer. Lung 198 (6), 897–907. doi:10.1007/s00408-020-00407-5
Araki, K., Youngblood, B., and Ahmed, R. (2013). Programmed Cell Death 1-directed Immunotherapy for Enhancing T-Cell Function. Cold Spring Harbor Symposia Quantitative Biol. 78, 239–247. doi:10.1101/sqb.2013.78.019869
Arner, E. N., Du, W., and Brekken, R. A. (2019). Behind the Wheel of Epithelial Plasticity in KRAS-Driven Cancers. Front. Oncol. 9, 1049. doi:10.3389/fonc.2019.01049
Božić, T., Kuo, C.-C., Hapala, J., Franzen, J., Eipel, M., Platzbecker, U., et al. (2021). Investigation of Measurable Residual Disease in Acute Myeloid Leukemia by DNA Methylation Patterns. Leukemia. doi:10.1038/s41375-021-01316-z
Cavalcante, G. M., Borges, D. P., de Oliveira, R. T. G., Furtado, C. L. M., Alves, A. P. N. N., Sousa, A. M., et al. (2020). Tissue Methylation and Demethylation Influence Translesion Synthesis DNA Polymerases (TLS) Contributing to the Genesis of Chromosomal Abnormalities in Myelodysplastic Syndrome. J. Clin. Pathol., jclinpath, 2020. doi:10.1136/jclinpath-2020-207131
Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017b). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019
Chen, Y. T., Shen, J. Y., Chen, D. P., Wu, C. F., Guo, R., Zhang, P. P., et al. (2020). Identification of Cross-Talk between m6A and 5mC Regulators Associated with Onco-Immunogenic Features and Prognosis across 33 Cancer Types. J. Hematol. Oncol. 13 (1), 22. doi:10.1186/s13045-020-00854-w
Chida, K., Kotani, D., Masuishi, T., Kawakami, T., Kawamoto, Y., Kato, K., et al. (2021). The Prognostic Impact of KRAS G12C Mutation in Patients with Metastatic Colorectal Cancer: A Multicenter Retrospective Observational Study. Oncol. 26 (10), 845–853. doi:10.1002/onco.13870
Cristall, K., Bidard, F.-C., Pierga, J.-Y., Rauh, M. J., Popova, T., Sebbag, C., et al. (2021). A DNA Methylation-Based Liquid Biopsy for Triple-Negative Breast Cancer. Npj Precis. Onc. 5 (1), 53. doi:10.1038/s41698-021-00198-9
Curran, W. J., Paulus, R., Langer, C. J., Komaki, R., Lee, J. S., Hauser, S., et al. (2011). Sequential vs Concurrent Chemoradiation for Stage III Non-small Cell Lung Cancer: Randomized Phase III Trial RTOG 9410. JNCI J. Natl. Cancer Inst. 103 (19), 1452–1460. doi:10.1093/jnci/djr325
Dor, Y., and Cedar, H. (2018). Principles of DNA Methylation and Their Implications for Biology and Medicine. The Lancet 392 (10149), 777–786. doi:10.1016/S0140-6736(18)31268-6
Du, J., Johnson, L. M., Jacobsen, S. E., and Patel, D. J. (2015). DNA Methylation Pathways and Their Crosstalk with Histone Methylation. Nat. Rev. Mol. Cel Biol 16 (9), 519–532. doi:10.1038/nrm4043
Fountzilas, E., Lampaki, S., Koliou, G.-A., Koumarianou, A., Levva, S., Vagionas, A., et al. (2021). Real-world Safety and Efficacy Data of Immunotherapy in Patients with Cancer and Autoimmune Disease: The Experience of the Hellenic Cooperative Oncology Group. Cancer Immunol. Immunother. doi:10.1007/s00262-021-02985-6
Gate, R. E., Cheng, C. S., Aiden, A. P., Siba, A., Tabaka, M., Lituiev, D., et al. (2018). Genetic Determinants of Co-accessible Chromatin Regions in Activated T Cells across Humans. Nat. Genet. 50 (8), 1140–1150. doi:10.1038/s41588-018-0156-2
Gautier, L., Cope, L., Bolstad, B. M., and Irizarry, R. A. (2004). Affy--analysis of Affymetrix GeneChip Data at the Probe Level. Bioinformatics 20 (3), 307–315. doi:10.1093/bioinformatics/btg405
Ghoneim, H. E., Fan, Y., Moustaki, A., Abdelsamed, H. A., Dash, P., Dogra, P., et al. (2017). De Novo epigenetic Programs Inhibit PD-1 Blockade-Mediated T Cell Rejuvenation. Cell 170 (1), 142–157. doi:10.1016/j.cell.2017.06.007
Grant, M. J., Herbst, R. S., and Goldberg, S. B. (2021). Selecting the Optimal Immunotherapy Regimen in Driver-Negative Metastatic NSCLC. Nat. Rev. Clin. Oncol. 18 (10), 625–644. doi:10.1038/s41571-021-00520-1
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Ito, S., Shen, L., Dai, Q., Wu, S. C., Collins, L. B., Swenberg, J. A., et al. (2011). Tet Proteins Can Convert 5-methylcytosine to 5-formylcytosine and 5-carboxylcytosine. Science 333 (6047), 1300–1303. doi:10.1126/science.1210597
Jazieh, A. R., Onal, H. C., Tan, D. S. W., Soo, R. A., Prabhash, K., Kumar, A., et al. (2021). Real-World Treatment Patterns and Clinical Outcomes in Patients with Stage III NSCLC: Results of KINDLE, a Multicountry Observational Study. J. Thorac. Oncol. 16 (10), 1733–1744. doi:10.1016/j.jtho.2021.05.003
Jiang, S. (2020). Tet2 at the Interface between Cancer and Immunity. Commun. Biol. 3 (1), 667. doi:10.1038/s42003-020-01391-5
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 28 (6), 882–883. doi:10.1093/bioinformatics/bts034
Li, B., Lu, Q., Song, Z. G., Yang, L., Jin, H., Li, Z. G., et al. (2013). Functional Analysis of DNA Methylation in Lung Cancer. Eur. Rev. Med. Pharmacol. Sci. 17 (9), 1191–1197.
Li, T., Pang, X., Wang, J., Wang, S., Guo, Y., He, N., et al. (2021). Exploration of the Tumor-Suppressive Immune Microenvironment by Integrated Analysis in EGFR-Mutant Lung Adenocarcinoma. Front. Oncol. 11, 591922. doi:10.3389/fonc.2021.591922
Lio, C. J., Yue, X., Lopez-Moyado, I. F., Tahiliani, M., Aravind, L., and Rao, A. (2020). TET Methylcytosine Oxidases: New Insights from a Decade of Research. J. Biosci. 45, 21. doi:10.1007/s12038-019-9973-4
Lussier, D. M., Alspach, E., Ward, J. P., Miceli, A. P., Runci, D., White, J. M., et al. (2021). Radiation-induced Neoantigens Broaden the Immunotherapeutic Window of Cancers with Low Mutational Loads. Proc. Natl. Acad. Sci. USA 118 (24), e2102611118. doi:10.1073/pnas.2102611118
Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501
Mayer, W., Niveleau, A., Walter, J., Fundele, R., and Haaf, T. (2000). Demethylation of the Zygotic Paternal Genome. Nature 403 (6769), 501–502. doi:10.1038/35000656
Miyakuni, K., Nishida, J., Koinuma, D., Nagae, G., Aburatani, H., Miyazono, K., et al. (2021). Genome‐wide Analysis of DNA Methylation Identifies the Apoptosis‐related Gene UQCRH as a Tumor Suppressor in Renal Cancer. Mol. Oncol. doi:10.1002/1878-0261.13040
Mo, Z., Cao, Z., Luo, S., Chen, Y., and Zhang, S. (2020). Novel Molecular Subtypes Associated with 5mC Methylation and Their Role in Hepatocellular Carcinoma Immunotherapy. Front. Mol. Biosci. 7, 562441. doi:10.3389/fmolb.2020.562441
Onodera, A., González-Avalos, E., Lio, C.-W. J., Georges, R. O., Bellacosa, A., Nakayama, T., et al. (2021). Roles of TET and TDG in DNA Demethylation in Proliferating and Non-proliferating Immune Cells. Genome Biol. 22 (1), 186. doi:10.1186/s13059-021-02384-1
Oswald, J., Engemann, S., Lane, N., Mayer, W., Olek, A., Fundele, R., et al. (2000). Active Demethylation of the Paternal Genome in the Mouse Zygote. Curr. Biol. 10 (8), 475–478. doi:10.1016/s0960-9822(00)00448-6
Pauken, K. E., Sammons, M. A., Odorizzi, P. M., Manne, S., Godec, J., Khan, O., et al. (2016). Epigenetic Stability of Exhausted T Cells Limits Durability of Reinvigoration by PD-1 Blockade. Science 354 (6316), 1160–1165. doi:10.1126/science.aaf2807
Prior, I. A., Lewis, P. D., and Mattos, C. (2012). A Comprehensive Survey of Ras Mutations in Cancer. Cancer Res. 72 (10), 2457–2467. doi:10.1158/0008-5472.CAN-11-2612
Qi, S., Li, Y., Dai, Z., Xiang, M., Wang, G., Wang, L., et al. (2019). Uhrf1-Mediated Tnf-α Gene Methylation Controls Proinflammatory Macrophages in Experimental Colitis Resembling Inflammatory Bowel Disease. J.I. 203 (11), 3045–3053. doi:10.4049/jimmunol.1900467
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Santha, S., Ling, X., Aljahdali, I. A. M., Rasam, S. S., Wang, X., Liao, J., et al. (2020). Mutant Kras as a Biomarker Plays a Favorable Role in FL118-Induced Apoptosis, Reactive Oxygen Species (ROS) Production and Modulation of Survivin, Mcl-1 and XIAP in Human Bladder Cancer. Cancers 12 (11), 3413. doi:10.3390/cancers12113413
Schübeler, D. (2015). Function and Information Content of DNA Methylation. Nature 517 (7534), 321–326. doi:10.1038/nature14192
Seremet, T., Koch, A., Jansen, Y., Schreuer, M., Wilgenhof, S., Del Marmol, V., et al. (2016). Molecular and Epigenetic Features of Melanomas and Tumor Immune Microenvironment Linked to Durable Remission to Ipilimumab-Based Immunotherapy in Metastatic Patients. J. Transl. Med. 14 (1), 232. doi:10.1186/s12967-016-0990-x
Siegel, R. L., Miller, K. D., and Jemal, A. (2020). Cancer Statistics, 2020. CA A. Cancer J. Clin. 70 (1), 7–30. doi:10.3322/caac.21590
Singhal, S. K., Byun, J. S., Park, S., Yan, T., Yancey, R., Caban, A., et al. (2021). Kaiso (ZBTB33) Subcellular Partitioning Functionally Links LC3A/B, the Tumor Microenvironment, and Breast Cancer Survival. Commun. Biol. 4 (1), 150. doi:10.1038/s42003-021-01651-y
Slieker, R. C., Roost, M. S., Van Iperen, L., Suchiman, H. E. D., Tobi, E. W., Carlotti, F., et al. (2015). DNA Methylation Landscapes of Human Fetal Development. Plos Genet. 11 (10), e1005583. doi:10.1371/journal.pgen.1005583
Smith, Z. D., and Meissner, A. (2013). DNA Methylation: Roles in Mammalian Development. Nat. Rev. Genet. 14 (3), 204–220. doi:10.1038/nrg3354
Treat, J. (2005). Incorporating Novel Agents with Gemcitabine-Based Treatment of NSCLC. Lung Cancer 50 (Suppl. 1), S8–S9. doi:10.1016/s0169-5002(05)81551-x
Wang, H., Chen, S., Meng, D., Wu, C., Zhu, J., Jiang, M., et al. (2021). Tumor Mutation burden and Differentially Mutated Genes Among Immune Phenotypes in Patients with Lung Adenocarcinoma. Ott. Vol. 14, 2953–2965. doi:10.2147/OTT.S294993
Weng, R. R., Lu, H. H., Lin, C. T., Fan, C. C., Lin, R. S., Huang, T. C., et al. (2021). Epigenetic Modulation of Immune Synaptic-Cytoskeletal Networks Potentiates γδ T Cell-Mediated Cytotoxicity in Lung Cancer. Nat. Commun. 12 (1), 2163. doi:10.1038/s41467-021-22433-4
Wu, C.-Y., Zhang, B., Kim, H., Anderson, S. K., Miller, J. S., and Cichocki, F. (2020). Ascorbic Acid Promotes KIR Demethylation during Early NK Cell Differentiation. J.I. 205 (6), 1513–1523. doi:10.4049/jimmunol.2000212
Wu, S. C., and Zhang, Y. (2010). Active DNA Demethylation: Many Roads lead to Rome. Nat. Rev. Mol. Cel Biol. 11 (9), 607–620. doi:10.1038/nrm2950
Wyatt, G. R. (1951). Recognition and Estimation of 5-methylcytosine in Nucleic Acids. Biochem. J. 48 (5), 581–584. doi:10.1042/bj0480581
Yoda, S., Dagogo-Jack, I., and Hata, A. N. (2019). Targeting Oncogenic Drivers in Lung Cancer: Recent Progress, Current Challenges and Future Opportunities. Pharmacol. Ther. 193, 20–30. doi:10.1016/j.pharmthera.2018.08.007
Yuan, M., Huang, L.-L., Chen, J.-H., Wu, J., and Xu, Q. (2019). The Emerging Treatment Landscape of Targeted Therapy in Non-small-cell Lung Cancer. Sig Transduct Target. Ther. 4, 61. doi:10.1038/s41392-019-0099-9
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020a). m6A Regulator-Mediated Methylation Modification Patterns and Tumor Microenvironment Infiltration Characterization in Gastric Cancer. Mol. Cancer 19 (1), 53. doi:10.1186/s12943-020-01170-0
Zhang, J., Zhong, X., Jiang, H., Jiang, H., Xie, T., Tian, Y., et al. (2020b). Comprehensive Characterization of the Tumor Microenvironment for Assessing Immunotherapy Outcome in Patients with Head and Neck Squamous Cell Carcinoma. aging 12 (22), 22509–22526. doi:10.18632/aging.103460
Keywords: lung adenocarcinoma, 5mC, tumour microenvironment, immunotherapy, mutation burden
Citation: Liu T, Guo L, Liu G, Hu X, Li X, Zhang J, Dai Z, Yu P, Jiang M, Wang J and Zhang J (2021) Molecular Characterization of the Clinical and Tumor Immune Microenvironment Signature of 5-methylcytosine-Related Regulators in non-small Cell Lung Cancer. Front. Cell Dev. Biol. 9:779367. doi: 10.3389/fcell.2021.779367
Received: 18 September 2021; Accepted: 19 October 2021;
Published: 11 November 2021.
Edited by:
Chunjie Jiang, University of Pennsylvania, United StatesReviewed by:
Weiru Liu, University of Pennsylvania, United StatesYang Xie, Brigham and Women’s Hospital, United States
Chong Jin, University of Pennsylvania, United States
Copyright © 2021 Liu, Guo, Liu, Hu, Li, Zhang, Dai, Yu, Jiang, Wang and Zhang. 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: Jian Zhang, zhangjian@gzhmu.edu.cn; Jian Wang, wangjian469@163.com
†These authors have contributed equally to this work