Molecular Characterization of the Clinical and Tumor Immune Microenvironment Signature of 5-methylcytosine-Related Regulators in non-small Cell Lung Cancer

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.
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.

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).

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 , 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 PCi 1 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: 5mCscore PCi 1 + PCi 2

Evaluation of Immune-Checkpoint Inhibitor Genomic and Clinical Information
To explore the application of the 5mC score to predict immunecheckpoint 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.

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).
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.

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 Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 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 The abundance of TME infiltrating cells in the 5mC modification patterns (*p < 0.05; **p < 0.01; ***p < 0.001; ns means not significant).
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 6 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).
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 7 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 immuneexcluded 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.

Identification of 5mC Methylation Gene Signature
To further identify the potential function of each m5C modification pattern, we determined 246 m5C phenotyperelated 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 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.
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 8 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 Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 9 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).

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 metacohort. 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 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.
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).
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.

DISCUSSION
DNA 5mC methylation is a dynamic and reversible posttranscriptional modification regulated by 5mC related regulators 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 phenotyperelated 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, DNMT3Amediated 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).
Frontiers in Cell and Developmental Biology | www.frontiersin.org November 2021 | Volume 9 | Article 779367 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 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.