Hyper-Methylated Hub Genes of T-Cell Receptor Signaling Predict a Poor Clinical Outcome in Lung Adenocarcinoma

Background Immune checkpoint inhibitors (ICIs) emerge as the first-line treatment of lung adenocarcinoma (LUAD); selection of subpopulations acquiring clinical benefit is required. Associations between epigenetic modulation of tumor microenvironment (TME) and clinical outcome are far from clear. We focused on immune-related genes closely regulated by DNA methylation to identify the potential clinical outcome indicators. Methods We systematically calculated immunophenotype score (IMpS) and classified immunophenotypes based on seven TME features in three independent cohorts. The overlapping of differential expressed genes and methylated probes targeted genes was regarded as genes closely regulated by DNA methylation. Then, probe/gene pairs which highly correlated with each other and IMpS were identified and named as immune-related probe/gene pairs (mIMg). Prognostic mIMg were selected and verified in seven independent validation cohorts. Results Three immune phenotypes were clustered, and similar results were obtained in the three independent training cohorts. C2 displayed as an immunologically hot phenotype, whereas C3 corresponded with immunologically cold phenotype. Average methylation level was decreased from C2 to C3 (C2 > C1 > C3). Similarly, ICIs nonresponders showed global hypo-methylation compared with responders. Genes in mIMg were mainly enriched, especially in T-cell receptor activation, and repressed in noninflamed TME by hyper-methylation. Among mIMg, low expression and hyper-methylation of CD247, LCK, and PSTPIP1 were risk factors of overall survival (OS). ICIs nonresponders were more likely to be hyper-methylated in the three genes. By integrating with the oncogenes status, we demonstrated that EGFR wt and SRGN overexpressed patients were associated with chronic inflammation and immune evasion, showing an immunologically hot phenotype, which might lead to the short OS but derive clinical benefit from ICIs. Conclusions This study identifies hyper-methylation and concurrent repression of CD247, LCK, PSTPIP1 as immune negative indicators and risk factors for prognosis in LUAD. Moreover, EGFR/SRGN axis may participate in immune modification to influence ICIs response and clinical outcome in LUAD.


Introduction
Nonsmall cell lung cancer (NSCLC) is the main subtype of lung cancer, which accounts for the primary cause of cancerassociated mortality. e treatment of lung adenocarcinoma (LUAD), the most common form of NSCLC, has achieved great progression in the last ten years with more and more targetable oncogene alternations identified [1]. What is more, immune checkpoint inhibitors (ICIs) targeting the programmed cell death protein 1 (PD-1) and its ligand (PD-L1) axis have markedly changed the first-line treatment of LUAD without molecular targets. Currently, extensive challenges still exist as only a subset of patients deriving clinical benefit [2].
Numerous studies had proved that the heterogeneity of tumor immune feature discriminates the clinical outcome of immunotherapies. e disturbance of immune effectors and regulators in the tumor microenvironment (TME) is just like the soil, fostering cancer cells (seeds) and participating in tumor development and progression. Dynamic equilibrium exists as the determinant of tumor immunogenicity, including the balance between infiltration of effectors (activated CD8+/CD4+ T-cells and Teffector memory (Tem) cells) and suppressors (regulatory T (Treg) cells and myeloid-derived suppressor cells (MDSCs)), the dominating type of T helper cells ( 1 and 2), and the expression of co-stimulators and co-inhibitors [3]. e immunophenotypes of TME predict response to ICIs. Immune-inflamed (also known as immunologically hot) TME was usually characterized by highly infiltrating T-cells accompanied by expression of suppressive molecules, such as PD-L1. e inflamed tumors profoundly respond well to ICIs. However, immune-excluded and immune-deserted (noninflamed or immunologically cold) TMEs are always associated with T-cell exclusion and unsusceptible to immunotherapy [4,5].
DNA methylation is a well-known epigenetic process, which participates in cancerogenesis and has been reported as a promising biomarker in the diagnosis of NSCLC [6,7]. With the widespread of ICIs, the relationship between TME and DNA methylation gradually grabbed much attention. TME shapes the fate of tumors with the aid of epigenetic alternations to modulate the dynamic gene expression. Moreover, DNA methylation silences or activates hub genes in the pivot immune signaling to shape TME and vice versa [8,9]. Although challenges exist, managing immune response to predict ICIs treatment outcome is a routine trend and outlook. So far, combinatorial analyses of epi-immune signatures have preliminarily declared the ability to predict survival and ICIs response [10]. Fewer hyper-methylated hub genes in critical ICIs response pathways correlated with improved ICIs response [11,12]. A simple prognostic profile or measure panel is an increasingly popular option in the routine use. However, even for the commonly and originally used biomarkers of ICIs response, such as PD-L1 status and tumor mutation burden (TMB), there are several exceptions for the differentiation between responders and non-responders. e predictive single hub genes in ICIs response merit further research. Hence, clinically available biomarkers for optimizing the use of ICIs and understanding the molecular determinants of immune response are still needed [13].
In this study, by integrating epigenetic and transcriptomic information of three LUAD cohorts, we explored the DNA methylation pattern of different immune phenotypes and identified hub genes of immune response which were regulated by DNA methylation and potentially influenced clinical outcome.
Gene length was calculated as the sum of lengths of nonredundant exon. Illumina nonnormalized summary-level data were read by "read.ilmn" function and normalized by the "neqc" function. Agilent single channel microarray intensity data were read by "read.maimages" function and processed by "backgroundCorrect" and "normal-izeBetweenArrays" unction. e abovedescribed processes were carried out by limma package [14]. e corresponding DNA methylation data (IDATs) including TCGA-LUAD, GSE56044, and GSE66836 methylation beta value were obtained [15,16]. Collectively, the three datasets were measured with Illumina Human-Methylation450 bead chip. Raw data were read by "read.metharray.exp" function to found rgSet [17] and filtered by "champ.filter" function. Specially, poor performing probes with detection p value more than 0.05, belonging to sex chromosome, known to have common SNPs at the CpG sites, or having been demonstrated to map to multiple places in the genome were removed prior to differential methylation analysis [18]. Normalization was performed by BMIQ method with "champ.norm" function. e DNA methylation pipeline was realized by minfi [17] and ChAMP [18] R package.
Five cohorts including GSE37745, GSE50081, GSE14814, GSE41271, and GSE42127 were utilized for the validation of prognosis associated genes [19][20][21][22][23]. Two cohorts with oncogene status including GSE13213, GSE11969 were used for verification analysis grouped by oncogene status [24,25]. DNA methylation cohorts containing ICIs response including GSE119144 and its corresponding expression data GSE135222 were used for the exploration of association between ICIs response and methylation/expression status. e raw files of Affymetrix were normalized by affy R package with Robust Multi-array Average (RMA) algorithm [26]. e procedures of processing Illumina and Agilent data were consistent with the training cohorts. Expression data from the same microarray manufactory were combined and the "ComBat" function in sva package was applied to adjust for batches deriving from different gene sets [27].
en, we weighted ssGSEA scores of positive immune factors (MHC molecules, effectors, 1 cells, co-stimulators) by "+1" and negative immune factors (suppressors, 2, co-inhibitors) by "−1" to calculate immune activity scores (IMaS) of every feature. Unsupervised hierarchical cluster was performed based on the IMaS of seven features and immunophenotype score (IMpS) was calculated by the sum of IMaS of every feature. e gene list containing 28 immune cell populations, MHC molecules, and immunostimulatory and immunoinhibitory factors was retrieved from the publication [3,30]. e consistency of different clusters was verified by cytolytic activity (CYT), which represented the ultimate effective mechanism in the cancer immunity cycle and was calculated as the geometric mean of granzyme A (GZMA) and perforin (PRF1) expression levels as previously defined [31]. e fraction of infiltrating cells were calculated using methylation data by EpiDISH [32]. e fraction of CD3+ T-cells (CD4+ T-cells and CD8+ T-cells) and epithelial cells were compared among different clusters.

Identification of Immune-Related Genes Closely
Regulated by DNA Methylation. Microarray probes annotation files were directly downloaded from GEO datasets. Homo. sapiens GRCh38.p13 GFF3 (v35) file downloaded from GENCODE website (https://www.gencodegenes.org/) was used for id transforming. Gene symbol was used as gene identification. If gene ID did not map one-to-one to the gene symbol, the first annotated gene was used to represent the others. RNA-seq data were preprocessed by normalizing distributions with "calcNormFactors" function in edgeR [33] and transforming to log2-counts per million (CPM) by the "voom" function in limma [14]. Linear modelling and empirical Bayes moderation in limma [14] were carried out to identify different expression genes (DEGs) of three datasets. Genes with pvalue less than 0.05 and the absolute value of log-2-fold-change (FC) more than 0 were considered as DEGs.
e average methylation level of each sample was calculated as mean of beta value of all methylation probes. Different methylation probes (DMPs) of three datasets were calculated using linear models implemented in ChAMP [18]. Probes with adjust pvalue less than 0.05 and the absolute value of log-2-fold-change (FC) more than 0 were considered as DMPs.
For assessing genes that underwent transcriptional regulation by gaining DNA methylation, we selected the overlapping of DEGs and DMPs corresponding genes, which were more likely to be different genes regulated by DNA methylation and named as mDEGs. Spearman correlation between methylation beta value and gene expression was carried out. Besides, the correlation between gene expression and DNA methylation beta value with IMpS was also performed. Probe/gene pairs whose correlation coefficient was no less than 0.5 and p value less than 0.05 were selected. e procedures described above were repeated in the three training cohorts, and the overlapping probe/genes were regarded as ultimate results. Probe/gene pairs which met the following requirements were regarded as the immune-related probe/gene pairs (mIMg): (1) the pairs were closely related with each other (|r| > 0.5, p.val <0.05); (2) both the probes and genes in the pairs were highly correlative with IMpS (|r| > 0.5, p.val <0.05. Beta value of each probe in mIMg was used for principal components analysis (PCA) in the three training cohorts by "PCA" function in FactoMineR package.

Functional Enrichment Analysis.
We further analyzed functions of the mDEGs. Gene Oncology (GO) functional enrichment was carried out in biological processes (BP), cellular components (CC), and molecular functions (MF). Terms with a pvalue < 0.05 were regarded as significant. e enrichment score was calculated as previously reported to combine the number and status (upregulated or downregulated) of genes enriched in the pathway [34]. Candidate genes were used for further stratification of overall population. Gene Set Enrichment Analysis (GSEA) was used to excavate mechanisms of them. Ontology gene sets in MSigDB were used for analysis [35]. Protein-protein interaction network (PPI) was performed by STRING website (https://string-db.org). Association between proteins encoded by genes in mIMg was represented by combined scores, which suggested the result of combinatorial analysis of gene neighborhood, fusion, occurrence, co-expression, and some other parameters [36]. R packages including clusterProfiler [37], GSVA, and ggplot2 were used for analysis and visualization.

Statistical Analysis.
Statistical analysis was conducted with R software (version 4, 4.0.4). e bioinformatic processes could be repeated with the R scripts uploaded on https://github.com/HU-ZX/01_mIMg_LUAD_2022. e datasets used in the study were showed in Table S1 and original data are accessible to download from the link attached in Github page. ANOVA test or Kruskal Wallis H test was used for continuous variables in multiple groups. T test or Mann Whitney U test was used for comparison between two groups. Differences between distributions of the groups were estimated by the Chi-squared test. Expression data and methylation beta value were directly used for the assessment of mIMg status. Overall, survival time (OS) and vital status of LUAD patients were used to represent clinical outcome. Patients were divided into two groups based on the threshold of each candidate marker. e thresholds were selected by "surv-cutpoint" function implemented in survminer R package to decrease the batch effect of calculation. Probes and genes in mIMg were, respectively, brought into multivariate Cox analysis by adjusting gender, age, stage to select potential prognostic biomarkers. e Kaplan-Meier method was used to estimate OS, and log-rank test between groups were performed. Hazard ratios (HRs) were derived from univariate Cox regressions of potential prognostic probe/ gene biomarkers. Statistical significance was set as p < 0.05. C-index based on univariate Cox model was calculated using the R package survcomp [38]. e receiver operating characteristic (ROC) curves at years 3, 5, and 9 years of the multivariate Cox model including age, gender, stage, and Journal of Oncology 3 candidate molecular factors were used to evaluate the discriminative ability of molecular indicators using the R package timeROC [39]. Efficiency of candidate molecular biomarkers was assessed by decision curve analyses [40].

Unsupervised Immunophenotype
Clusters. e workflow of this study is shown in Figure 1. Basic information and clinical pathological features of training and validation cohorts used in this study are showed in supplementary tables Table S2 and Table S3, respectively. ere were 450 patients in the TCGA cohort, 112 patients in the GSE66863 & GSE66836 cohort, 78 patients in the GSE60644 & GSE56044 cohort for the integrative analysis of transcriptomic and epigenetic data. Samples were clustered into three clusters based on the IMaS of 7 immune features respectively in the above-described cohorts. Similar cluster results were obtained in different datasets and showed by the heatmaps (Figures 2(a), 2(e), and 2(i)).
Both immune effectors and regulators were in the dominant position in the C2 cluster. In detail, tumors of the C2 cluster were characterized by the presence of CD8+ and CD4+ T-cells which act as defenders in the proximity of tumors. However, the immunosurveillance system in C2 was also activated which manifested as high expression of co-inhibitors including immune checkpoints and other immunosuppressed factors. C2 cluster presented as preexisting anti-tumor immune response, tending to be immunologically hot phenotype, and was more likely to benefit from immunotherapies.
In contrast, the C3 cluster was characterized by fewer infiltrated TILs and low expression of immune regulators, which was more likely to be immunologically cold or noninflamed tumor. Not surprisingly, tumors in the C3 cluster belonged to ICIs-resistant subset. C1 cluster demonstrated high immune heterogeneity and could not be classified into C2 or C3 cluster (Figures 2(a), 2(e), 2(i)).

Methylation Features of Different Immunophenotypes.
e average beta value of each sample was used for assessing the global methylation features. Notably, the average methylation level was consistent with the trend of IMpS (C2 > C1 > C3, Kruskal-Wallis test, p � 2.7e − 12, p � 0.00028, p � 4.1e − 05, Figures 3(a)-3(c)), which validated the result published before [9]. e same trend was observed in the ICIs responders and nonresponders, albeit no significant difference was observed (Figure 3(d)).
C3 (noninflamed tumor cluster) and C2 (inflamed tumor cluster) had different immune features; 4847 overlapping probe/gene pairs of DEGs and DMPs among the three cohorts were identified (Figures 3(e)-3(h)). 3806 probes were hypo-methylated in a noninflamed tumor. e cisregulatory pattern was defined as hypo-methylation repressed gene expression or hyper-methylation activated gene expression, in which condition, the changing trend of gene expression and DNA methylation level was the same. e trans-regulatory pattern was defined as hyper-methylation repressed gene expression or hypo-methylation activated gene expression. In our study, we confirmed that DNA methylation was more likely to trans-regulate gene expression (Chisq test, X-squared = 25, p � 2.035e − 07, Figure 3(i)). DNA methylation was unevenly distributed. TSS1500 contained more hypo-methylated probes in C3 compared with non-TSS1500 region (Chisq test, X-squared = 28.08, p � 1.164e − 07, Figure 3(j)). Probes located in promoter (TSS200, TSS1500, 1 st exon, 5'UTR) were more likely to trans-regulate gene expression while those located in gene body tended to cis-regulate gene expression (Chisq test, X-squared = 202.34, p < 2.2e − 16, Figure 3(j)). Flanking regions of CG island (CGI), also known as "shore," were enriched with more trans-regulative probes compared with CGI region (Chisq test, X-squared = 51.129, p � 8.648e − 13, Figure 3(k)). e chromosomal distribution characteristics of probes were presented in Figure 3(l).

Hyper-Methylation of Hub Genes in Immune Response
Predicted Noninflamed Phenotype. Functional enrichment analysis showed that the mDEGs enriched in T-cell activation and some other immune-related processes. Besides, it suggested GTPase and Src homology 2 (SH2) domain as potential targets of DNA methylation regulatory system in immunophenotype. GTPase is a rate-limiting enzyme in GTP hydrolysis to participate in G protein inactivation. SH2 domain is a special structure in the combination of enzymes and prolongs docking time. e two are in the downstream of TCR activation signaling. e processes described above were all downregulated in the C3 cluster (Figure 4(a)).
Spearman correlation identified 168 probes and 86 corresponding genes which highly correlated with IMpS (Supplementary Table S4). Meanwhile, correlation between methylation level and gene expression was performed, and 61 probes whose beta level were highly correlated with its corresponding gene expression were selected (Supplementary Table S5). Ultimately, 24 probes and 19 corresponding genes which were contained in both two probe/gene sets described above were defined as mIMg (Supplementary  Table S6, Figures 4(b) and 4(c)). PCA analysis of three training cohorts indicated that mIMg successfully distinguished patients with different immunophenotypes in the three LUAD training cohorts (Figures S1(a)-S1(c). Trans-regulative pattern was more common in mIMg and genes in mIMg tended to be repressed by hyper-methylation in the noninflamed tumor (Supplementary Table S6 Methylation level of IMpS-related mDEGs according to ICIs response in the GSE119144 cohort were investigated, and 44 of 168 probes were significantly different between responders and nonresponders with the same trend in the training cohort (Supplementary Table S7-1). 10 probes were in mIMg, including cg07786657, cg09032544, cg07728874, cg24841244, cg11384427, cg08450017/, cg11683242, cg14145194, cg15518883, and cg26227523. However, expression of the corresponding genes was not significantly different (Supplementary Table S7-2).

Silence of T-Cell Activation Associated Genes by DNA Methylation Predicted a Poor Clinical Outcome.
To investigate the prognostic capacity of mIMg, multivariate Cox analyses adjusted for age, gender, stage based on hyper/hypo methylation (high/low expression) groups of the probe/gene pairs were performed, respectively, in the TCGA-LUAD cohort, and 15 probes together with their corresponding 13 genes distinguished OS (Supplementary Table S8). Prognostic ability of genes in mIMg was verified in two integrative expression cohorts (Supplementary Table S3). ArfGAP with coiled coil, ankyrin repeat and PH domains 1 (ACAP1), Rho GTPaseactivating protein 30 (ARHGAP30), CD247, lymphocyte transmembrane adaptor 1 (LAX1), lymphocyte-specific protein tyrosine kinase (LCK), proline-serine-threonine phosphatase interacting protein 1 (PSTPIP1) stood the test (Table 1). Among these prognostic mIMg factors, hypermethylation and concurrent downregulation of gene expression was significantly blamed for low OS.

Crosstalk between DNA Methylation and Oncogenic
Mutations. Growing evidence suggests that Epidermal Growth Factor Receptor (EGFR) and other oncogenic driver mutations modify the TME with low immune checkpoints and TMB thus reduce ICIs response [41]. More mechanisms associated with ICIs resistance are gradually surfacing and DNA methylation is one. Average methylation level of three clusters stratified by EGFR status, KRAS status, and tumor protein p53 (TP53) status was consistent   Multivariate Cox analyses were carried out according to different oncogene status, respectively, in TCGA-LUAD cohort. IMpS-related mDEGs whose expression and methylation level concurrently had prognostic ability were used for further multivariate Cox validations in GSE11969 and GSE13213 cohorts (Supplementary Table S9). Although prognostic mIMg was found in EGFR mut, KRAS mut/wt, TP53 mut/wt group, these genes did not stand validation tests. mDEGs which potentially influenced prognosis were presented in the forest plots (Supplementary Table S9, Figures S7(G)-S7(I)). cg02851793/Serglycin (SRGN) was a prognostic pair according to EGFR status. Expression of SRGN was regulated by methylation regardless of EGFR status (Figures 7(a) and 7(b)). e total group was stratified into SRGN high and SRGN low groups based on OS. Low expression of SRGN and hyper-methylation in cg02851793 was a protective factor in the EGFR wt group of TCGA-LUAD cohort (Figures 7(c)-7(e)). Low expression of SRGN also predicted a better OS in the EGFR wt group of GSE11969 and GSE13213 cohorts (Figures 7(c), 7(f ) and 7(g)). However, low expression of SRGN was a risk factor in the EGFR mut group in TCGA-LUAD and GSE11969 but failed to be validated in GSE13213 (Figure 7(c)). Methylation level of SRGN did not affect survival in the EGFR mut group. Log-rank test was used to compare OS between SRGN overexpressed and low-expressed group. C-index based on univariate model of SRGN was calculated in EGFR wt patients (Figures 7(d)-7(g)). ROC curves based on multivariate cox model were used to evaluate the prognostic ability in training and validation cohorts (Figures 7(h)-7(k)). We demonstrated that the expression level of SRGN cooperating with EGFR status could discriminate clinical outcome. en, we explored the expression of SRGN and EGFR status. Although comparison of SRGN in EGFR mut and wt groups was not significant in TCGA cohort and GSE13213 validation cohort (p � 0.7341, p � 0.09173), SRGN was higher expressed in EGFR wt group in both SRGN high and SRGN low groups (Figures 7(l), 7(n)). e expression SRGN and EGFR tended to be low correlated with each other, yet it was not validated successfully (Figures 7(m), 7(o)). GSEA analysis was carried out to identify potential mechanisms. Notably, all patients in EGFR wt and SRGN high group shared inflamed phenotypes (C1 and C2). Immune effective cells were highly infiltrating accompanied by expression of immune inhibitors in this group (Figures 8(a) and 8(b)). MHC-I antigen presentation process and 1 cells activity were enriched in this group (Figure 8(d)). Moreover, initiate immune pathways were also enriched, indicating chronic inflammation (Figure 8(d)). In contrast, EGFR mut and SRGN low group was characterized as a noninflamed phenotype. EGFR wt and SRGN low group showed immune heterogeneity and was well classified by immune cluster. In addition, SRGN was higher in ICIs responders compared with nonresponders, which confirmed that SRGN might function to modify TME potentially (Figure 8(c)). To sum up, with the involvement of immune disturbance, EGFR wt and SRGN high patients were associated with short OS but might benefit from ICIs.

Discussion
In recent years, ICIs have achieved impressive success in the treatment of LUAD. Selection of responders is becoming a more and more crucial issue. ere is a growing awareness that cancer cells are fostered by highly heterogeneous and plastic cells in TME engaging in well-orchestrated reciprocal interactions [42]. Immunologically hot or inflamed TME is characterized with high density of immune effectors accompanied by the expression of immune checkpoints. In this case, ICIs reinvigorate anti-tumor immune response. In contrast, immunologically cold or noninflamed TME lacks ICIs response [4]. In our study, three immune phenotypes were classified according to seven immune features; C2 (k)  cluster corresponded to immunological hot phenotype, indicating preexisting anti-tumor immune response, and was more likely to benefit from immunotherapies. C3 cluster had fewer infiltrated TILs accompanied by low expression of immune regulators, tending to be an immunologically cold or noninflamed tumor. Tumors in C3 cluster were more likely to be ICIs-resistant subset. DNA methylation is a key regulatory strategy in epigenetic modulation of gene expression. It was previously reported that hyper-methylation on the promoters of tumor suppressor genes could promote tumorigenesis [43,44]. Emerging evidences have suggested the critical role of DNA methylation in immune modification of tumors and mediating the immune response [8,9] . Recent studies had demonstrated that DNA methylation loss counteracted with TMB and copy number load formed by genome instability in mitotic cell division and involved in immune evasion to increase ICIs resistance [45]. We confirmed that noninflamed TME and ICIs nonresponders underwent methylation loss. e result implies that the combinatorial regimen of DNA methyl-transferase inhibitors (DMTis) and ICIs holds promise for improving the therapeutic efficacy of ICIs [46].
At present, biomarkers are continuing to be excavated for directing ICIs response as target therapy. Emerging studies are developing prediction models for immune response or clinical outcome [4]. Herein, we focused on genes involving in immune modification which are closely regulated by DNA methylation, naming these methylation probe/gene pairs as mIMg. Although the average methylation level was low in the immunologically cold cluster and ICIs nonresponders, hub genes in immune-related pathways, especially in T-cell receptor activation, were silenced by hyper-methylation.
As we all know, two distinct steps must be completed to elicit an effective anti-tumor response. Firstly, dendritic cells (DCs) accomplish antigen cross-presentation for CD8+ T-cell initiation. en, antigens must be directly presented by the tumor via MHC class I pathway and recognized by activated CD8+ T-cells for killing [5]. Among mIMg, most genes in noninflamed tumor were repressed by hypermethylation and involved in CD8+ T-cell activation, especially in the process of interaction between MHC class I and T-cell receptor (TCR) complex. TAP1 encodes antigen peptide transporter 1, which involves in the transport of antigens from the cytoplasm to the endoplasmic reticulum for association with MHC class I molecules. CD247, CD3D encode zeta and delta chains of CD3 which participate in TCR complex. Zeta chains transduce positive signals by immunoreceptor tyrosine-based activation motifs (ITAMs). Apart from editing of TCR, the downstream signal also experienced hyper-methylation, such as LCK, which encodes Syk family kinases Lck and associates with CD4 or CD8 cytoplasmic tails. LCK phosphorylates the tyrosine residues of ITAMs in CD3 and zeta chains and transduces the TCR activation signal in an SH2 domain-dependent manner [47]. Moreover, the hyper-methylation of costimulators and regulators in the CD8+ T-cell activation such as intercellular adhesion molecule 3 (ICAM3), PSTPIP1, CD27, CD37, signaling threshold regulating transmembrane adaptor 1 (STI1), LAX1, Myosin IG (MYO1G) also proved to participate in immunologically cold phenotype. We then compared methylation/expression level of mIMg according to ICIs response for validation, and found that the methylation level of CD247, CD3D, LCK, ICAM3, STI1, PSTPIP1, and CXCR6 were enhanced in the nonresponders, which verified our discovery. Even so, no significant difference was found in the expression of genes described above. As 58 samples in GSE119144 had complete methylation data whereas only 27 samples in the corresponding expression cohort GSE135222 had transcriptomic data, we speculated that the failure of distinction between  responders and nonresponders in transcriptomic level could come down to the small sample size. Moreover, as epigenetic modulation occurred earlier than expression change, methylation changes might be more sensitive than expression changes.
To investigate the clinical value in genes repressed by hyper-methylation in noninflamed TME, multivariate Cox analyses adjusted for universally known prognostic factors (gender, age, stage) were carried out among mIMg. We demonstrated that repression of CD247, LCK, PSTPIP1 by hyper-methylation corresponded with immunological cold phenotype and were risk factors of OS. e crucial role of CD247 and LCK need not to be restated more. PSTPIP1 encodes a CD2 cytoplasmic tail-binding protein and acts as an immunoinhibitor by blocking CD2/CD58 contact. CD2 interacts with CD58 at the initial stages of T-cell activation even prior to the recognition of TCR and MHC by increasing the dwell time between antigen presentation cells (APC) and T-cell [48]. In our study, we suggested that hyper-methylation and concurrent repression of CD247, LCK, and PSTPIP1 attenuated CD8+ T-cell activation by disturbing MHC I antigen binding, TCR signal transduction, and  regulation of co-stimulators. In this way, noninflamed TME formed with inferior ICIs response and bad OS. Apart from hyper-methylation of hub genes in T-cell activation, we also identified that hyper-methylation and downregulation of ARHGAP30 and ACAP1 as potential risk elements related to noninflamed TME. Although both methylation and expression level had no significant difference according to ICIs response, the important role of them in immunophenotype distinction could not be denied. e two elements both encode GAPs, which act by binding to the GTPase and promoting the conversion of small guanine nucleotide-binding proteins (small G proteins) from active GTP bound state to inactive GDP bound state. e representative substrate of GTPase are Ras family proteins, which are activated for participating in proliferation, attachment, motility, and promote malignancy. ARHGAP30 and ACAP1 are both GAPs for Ras superfamily proteins and past studies have demonstrated that over-expression of them attenuated malignant characteristics and acted as favorable prognostic factors [49,50]. In our study, we obtained similar results and suggested that over-expression of the two genes were closely related to DNA hypo-methylation. As GAPs are downstream signal of T-cell activation, we supported the opinion that hypo-methylation of ARHGAP30 and ACAP1 were the outcomes of inflamed TME, instead of the reason.
In most of the ICIs clinical studies, patients that harbored oncogenic alternations were excluded [39]. Past studies showed that patients with sensitive mutations had lower TMB and immunogenicity compared with those wild-type patients [51]. Coinciding with the trend of the overall population, we found that the difference of the average methylation level was more pronounced according to immunophenotypes in the EGFR wt and KRAS wt subgroups. However, a smaller sample size of patients with oncogene mutation might also make sense. Methylation features of patients with or without TP53 alternations were similar. TP53 is a pro-oncogene, whose alternations are in the prophase of tumorigenesis and earlier than other epigenetic adaptation changes. In our study, we recognized EGFR/SRGN axis as a potential mechanism to discriminate ICIs response and OS. SRGN encodes a kind of proteoglycan, which can be secreted by the extracellular matrix of tumor cells to create a pro-inflammatory TME and is regarded as a driving factor of aggressive phenotype [52]. In our study, we illustrated that the over-expression of SRGN was closely regulated by hypo-methylation despite the EGFR status but the EGFR wt group tended to express more SRGN. EGFR wt patients accompanied with high expression of SRGN correlated with short OS but might benefit from ICIs. We demonstrated that EGFR wt with SRGN over-expression patients was displayed as inflamed tumors, but associated with "cancer-promoting inflammation", which shaped TME toward a tumor-permissive state by chronic inflammation and immune evasion [53].
Our study had some limitations. Firstly, although we declared that mIMg played a crucial role in the immune response, and CD247, LCK, PSTPIP1, and SRGN were potential prognostic indicators, further analyses and largescale clinical validations are needed. e association between mIMg and ICIs response deserves far more confirmations.
Secondly, the measure of methylation beta value and expression data of every cohort might be influenced by many factors; in addition to samples' attributes and testing conditions, the threshold of methylation and expression level of different genes had distinct effects on results. In this condition, we set the cutoff value according to OS in every cohort, respectively, to reduce the batch effect. A consistent and steady cutoff value in the overall population may be needed. Besides, although prognostic ICIs response-related indicators were recognized, the direction of the association between immune response and methylation/ expression status could not be determined. Past study demonstrated that ICIs could reshape the TCR repertoire, reiterating that TME was plastic and in continuous modification [54]. Hyper-methylation of CD247, LCK, and PSTPIP1 might be a way of TCR editing to influence immune response and survival. However, we could not determine whether it was hypermethylation of TCR genes that led to immune exclusion or the immune surveillance that shaped or incapacitated TCR by DNA methylation. Moreover, oncogene mutations are numerous and in the dynamic change. We did not bring oncogenic alternations into multivariate Cox analysis in the process of choosing prognostic variables. To eliminate the effect of multicollinearity among the explanatory variables, we performed an extra hierarchical analysis according to oncogene status.

Conclusions
is study demonstrated the pivotal role of DNA methylation in immune response in LUAD. TCR editing by hypermethylation of CD247, LCK, and PSTPIP1 acts as potential immune response indicators and prognostic factors. EGFR/ SRGN axis involves in TME modification to influence clinical outcomes.

Data Availability
Publicly available datasets were analyzed in this study. e details and downloading websites of these data can be found in Table S1. e bioinformatic process is accessible on Github (https://github.com/HU-ZX/01_mIMg_LUAD_2022).

Ethical Approval
e study was carried out with the ethical approval of the institutional Ethics Committee of the Faculty of Medicine at China-Japan Friendship Hospital approval (2018-99-K71). Publicly available data were used in the article, so the ethical approval was waived.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
ZX-H and HJ-C contributed to the idea and design. XY-L, J-L, X-Z, YX-Y, and KX-T hunted for and downloaded the data. ZX-H, CX-X, and HJ-D contributed to data excavation. ZX-H, HJ-C, and JB-Z accomplished the manuscript writing and revision. All authors have read and approved the submitted version.

Supplementary Materials
Table S1: Information of data was utilized in our study. Table  S2: Clinical pathological features of the training cohorts. Table  S3: Clinical pathological features of validation cohorts. Table  S4: Methylation probes and corresponding regulating genes which both are highly correlated with IMpS. Table S5: Methylation probes/gene pairs which are closely correlated with each other. Table S6: Features of mIMg in immune response in the training cohorts. Table S7-1: Comparison of the methylation level of probes highly correlated with IMpS between ICIs responders and nonresponders. Table S7-2: Comparison of the expression level of corresponding genes in Table S7-1 between ICIs responders and nonresponders. Table S8: Multivariate Cox analysis of methylation probes in mIMg pairs which both distinguished OS in the TCGA-LUAD cohort. Table S9: Prognostic methylation probes according to oncogenes status. Figure S1: PCA analysis showed the distinct methylation status of genes in mIMg in three immunophenotypes. mIMg was defined as DNA methylation probe/gene pairs, genes in which were immune response-related hub genes and closely regulated by DNA methylation. Figure S2: Correlation between cg09032544/ CD247, cg07786657/CD247, cg11683242/LCK, cg26227523/ PSTPIP1 in mIMg and 7 immune features in GSE60644 & GSE56044, GSE66863 & GSE66836 cohorts. Figure S3: Evaluation of CD247, LCK, and PSTPIP1 in predicting ICIs response. (A-C) No significant difference was observed in the comparison of CD247, LCK, and PSTPIP1 expression levels between ICIs responders and nonresponders in GSE135222. Figure