5-methylcytosine RNA methylation regulators affect prognosis and tumor microenvironment in lung adenocarcinoma.

Background: Accumulating evidence has shown that 5-methylcytosinec (m5C) RNA methylation plays an essential role in tumorigenesis. However, the roles of m5C regulators in the prognosis, tumor microenvironment (TME), and immunotherapy responses of lung adenocarcinoma (LUAD) have not been fully analyzed.
Methods: Based on 14 m5C RNA regulators, we evaluated the m5C RNA modification patterns in patients with LUAD (n=594) in The Cancer Genome Atlas (TCGA). Unsupervised clustering analysis was performed to confirm distinct m5C modification patterns. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to investigate the biological functions of differentially expressed genes (DEGs) among different m5C RNA modification patterns. An m5C signature (m5Csig) was constructed using least absolute shrinkage and selection operator (LASSO) algorithms. The GSE72094 cohort (n=442) from the Gene Expression Omnibus (GEO) was used to validate m5Csig. A receiver operating characteristic (ROC) model was constructed to evaluate the sensitivity and specificity of m5Csig. Tumor-infiltrating immune cells (TIICs) between the high- and low-risk groups were estimated using the Cell Type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm.
Results: We identified 3 m5C RNA modification clusters. Overall survival (OS) differed among the 3 clusters. The m5Csig, including TRDMT1, NSUN1, NSUN4, NSUN7, and ALYREF, was constructed to classify patients with LUAD into high- and low-risk groups. The high-risk group, with more immune cell infiltration, had a significantly poorer OS than that the low-risk group, which was associated with better response to immune checkpoint blockade therapy. Conclusions: The present study revealed that m5C RNA regulators play a significant role in TME regulation in LUAD. The m5Csig can predict the prognosis of patients with LUAD and might provide novel strategies for tumor immunotherapy. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer is one of the leading causes of cancer-related morbidity and mortality worldwide (1). Non-small cell lung carcinoma (NSCLC), accounting for 85% of all lung cancers, mainly comprises lung squamous cell carcinoma (LUSC) and lung adenocarcinoma (LUAD), with LUAD being the most common NSCLC subtype (2). With advances in targeted therapy and immunotherapy (3,4), the 5-year survival rate of patients with NSCLC is still unsatisfactory, at 4–17% (5). Patients with the same clinical characteristics can have distinctly different prognoses because of molecular differences. Therefore, there is an urgent need to confirm new molecular targets to improve the clinical treatment outcome in patients with LUAD.

Methylation of RNA, an important epigenetic modification that includes 5-methylcytosine (m5C), N6-methyladenosine (m6A), N1-methyladenosine (m1A), pseudouridine (Ψ), and inosine (I), has been identified to decorate protein-coding messenger RNAs (mRNAs) and noncoding RNAs (ncRNAs) (6-10). Modifications of RNA play crucial roles in RNA translation, transcription, processing, stability, and splicing (11,12), and m5C is one of the most common RNA modifications (13). The m5C RNA methylation can be catalyzed dynamically by a series of significant mediator proteins known as “writers” [tRNA aspartic acid methyltransferase 1 (TRDMT1), NOP2 nucleolar protein (NSUN1), NOP2/Sun RNA methyltransferase 2 (NSUN2), NSUN3-7], “readers” [Aly/REF export factor (ALYREF) and Y-box binding protein 1 (YBX1)], and “erasers” [tetmethylcytosine dioxygenase 1 (TET1), TET2-3]” (13-16). Dysregulation and disorder of m5C are associated with the occurrence of human diseases, including malignancies (17-19).

In recent years, immune checkpoint blockade (ICB) has made great breakthroughs in clinical efficacy for patients with cancer (20). However, only a small number of patients benefit from ICB (21). Numerous studies have identified that the tumor microenvironment (TME), which contains immune cells (such as T and B lymphocytes, natural killer (NK) cells, macrophages, polymorphonuclear cells, dendritic cells, as well as mast cells) and stromal cells, plays a crucial role in tumor progression, immunotherapy response, and immune escape (22,23). However, the relationship among m5C, immunotherapy response, and the TME in LUAD remains unclear. Therefore, a comprehensive understanding of the effect of m5C regulators on the TME might provide new insights into the immune regulation of the TME.

In this study, we analyzed 14 m5C RNA methylation regulators in LUAD from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. We identified that m5C regulators were closely associated with LUAD prognosis, and then constructed an m5C signature (m5Csig) to predict the LUAD survival and evaluate the response to ICB. Overall, the results indicated that m5Csig could act as a biomarker to predict survival and the response of ICB in LUAD. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-500/rc).

Methods

Acquisition of data

The RNA sequencing data (n=594), somatic mutation information (n=561), copy number variation (CNV) information (n=555), and the corresponding clinicopathological features (n=522) of patients with LUAD were downloaded from TCGA (https://portal.gdc.cancer.gov/). To validate the findings in the TCGA database, the validation cohort GSE72094 (n=442) was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Protein-protein interaction analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) was used to analyze protein-protein interaction (PPI) information and detect the interactions of 14 m5C regulators. We then extracted PPI pairs with a combined score of 0.4.

Unsupervised clustering for 14 m5C regulators

Unsupervised clustering analysis was used to identify different m5C RNA modification patterns among patients with LUAD (n=535) from TCGA. The 14 m5C regulators included 8 writers (TRDMT1, NSUN1–7), 2 readers (ALYREF and YBX1), and 4 erasers [TET1–3, AlkB homolog 1, histone H2A dioxygenase (ALKBH1)]. The consensus clustering algorithm was employed to categorize patients with LUAD into different modification patterns (24). The consensus ClusterPlus package (https://bioconductor.org/packages/release/bioc/html/ConsensusClusterPlus.html) was used to perform the above steps with a cycle computation of 1,000 iterations to guarantee the stability and reliability of the results (25). The overall survival (OS) rates of patients with the 3 modification patterns were calculated using the Kaplan-Meier method.

Identification of differentially expressed genes between m5C modification patterns

The empirical Bayesian approach in the limma R package (https://bioconductor.org/packages/release/bioc/html/limma.html) was applied to identify differentially expressed genes (DEGs) among the different m5C modification patterns in the standard comparison mode (26). The significance criteria for determining DEGs was set as an adjusted P value <0.001. To identify the potential functions and pathways enriched in the different modification patterns, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed based on the DEGs identified among the different modification patterns (27).

Development of the m5C regulators-related prognostic signature

Univariate Cox regression analysis was performed to identify m5C RNA methylation regulators that were associated with the OS of patients with LUAD. A least absolute shrinkage and selection operator (LASSO) Cox regression algorithm was carried out to construct an optimal m5Csig to predict the prognosis of patients with LUAD. Then, the obtained prognosis-associated genes were used to construct the m5Cscore function to calculate the score for each patient. The m5Cscore formula we used as follows:

where m5Cscore is a prognostic risk score for patients with LUAD patients, coefi represents the coefficient, and expri represents the expression of each prognostic gene. According to the median m5Cscore, the patients with LUAD were classified into high- and low-risk groups. The Kaplan-Meier method with the log-rank test was used to evaluate the OS differences between the high- and low-risk groups, and a receiver operating characteristic curve (ROC) was used to evaluate the prediction accuracy of m5Csig. Univariate and multivariate Cox regression analyses were used to explore prognostic values of m5Csig and clinical characteristics.

Estimation of TME cell infiltration

According to the method used by Newman et al. (28), 22 types of tumor-infiltrating immune cells (TIICs) between the high- and low-risk groups were estimated using the Cell Type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm.

Building and validation of a nomogram

Clinical factors [age, gender, stage and tumor-node-metastasis (TNM) stage] and the m5Cscore were used to develop a prognostic nomogram to evaluate the probability of 1-, 2-, and 3-year OS for patients with LUAD (29). The C-index and a calibration plot were constructed to estimate the accuracy and consistency of the m5Cscore.

Statistical analysis

Spearman and distance correlation analyses were used to compute the correlation coefficients among the expression levels of m5C regulators. The Wilcoxon test was used to analyze the difference between 2 groups, and the Kruskal–Wallis test and one-way analysis of variance (ANOVA) were used among 3 or more groups. LASSO Cox regression and Kaplan-Meier analyses were performed to construct and evaluate the m5Cscore. The area under the receiver operating characteristic curve (AUC) was used to investigate the time-dependent prognostic value of the m5Cscore. Multivariate Cox regression and stratified analysis were used to verify the independence of the m5Cscore from other clinical factors. Statistical significance was set at P<0.05, and all statistical P values were 2-sided. All data were processed using R 4.0.3 software (R Foundation for Statistical Computing, Vienna, Austria).

Results

Landscape of genetic variation among m5C regulators in LUAD

The workflow of our study is shown in Figure S1. A total of 14 m5C regulators including 8 writers (TRDMT1, NSUN1–7), 2 readers (ALYREF and YBX1), and 4 erasers (TET1–3, ALKBH1) were included in our study (Table S1). First, the frequency of CNVs and somatic mutations of the 14 m5C regulators were investigated in LUAD. The CNV alteration frequency indicated that CNV alterations were ubiquitous among the 14 m5C regulators. As shown in and Table S2, NSUN2 (13.69% amplification vs. 1.80% deletion), ALYREF (10.81% amplification vs. 1.44% deletion), YBX1 (7.39% amplification vs. 1.80% deletion), and NSUN4 (6.13% amplification vs. 2.52% deletion) were associated with amplification of the copy number, while ALYBH1 (4.86% deletion vs. 1.80% amplification) and NSUN1 (5.95% deletion vs. 4.32% amplification) were frequently deleted. The distribution of CNV alteration of m5C regulators on chromosomes is shown in . The analysis showed that 13.19% of patients with LUAD (n=74) experienced mutations of m5C regulators. The highest mutation frequency was exhibited by TET1 (4%) followed by TET2 (2%) and TET3 (2%), while the genes including the writers (NSUN1, NSUN3–5 and NSUN7), readers (ALYREF and YBX1) had no mutations in the patients with LUAD (). To explore the relationship between CNV alteration and the expression of m5C regulators, we analyzed the mRNA expression levels of the regulators. The results indicated that the expression levels of NSUN1, NSUN2, NSUN4–7, ALKBH1, TET1, TET3, and ALYREF were significantly upregulated in LUAD (P<0.001), whereas the expression level of TRDMT1 was significantly downregulated in LUAD (P<0.001). No significant difference was found in the expression levels of TET2, YBX1, and NSUN3 (). These analyses showed CNV might play a crucial role in the imbalanced expression of m5C regulators, which could affect the occurrence and progression of LUAD. The clinicopathological features of the patients with LUAD are summarized in Table S3.

Figure 1

Landscape of the genetics and expression analysis of m5C regulators in TCGA-LUAD. (A) The CNV frequency of 14 m5C regulators; (B) the distribution of CNV alterations of the 14 m5C regulators on 23 chromosomes; (C) the mutation frequency of the 14 m5C regulators in 561 patients with LUAD; (D) the expression levels of the 14 m5C regulators between LUAD tissues and normal tissues; (E) heatmap of the 14 m5C regulators between LUAD tissues and normal tissues; (F) correlations among the 14 m5C regulators; (G) PPI network of the 14 m5C regulators; (H) the interaction numbers of each regulator with the other 13 regulators. ***P<0.001. TCGA-LUAD, The Cancer Genome Atlas lung adenocarcinoma cohort; CNV, copy number variation; PPI, protein-protein interaction.

Correlation and interaction between m5C regulators

To determine the crosstalk among m5C regulators, correlation analysis was performed, which showed mainly positive correlations among the m5C regulators; however, several regulators exhibited a negative correlation among the m5C regulators. For example, TET2 and NSUN3 had the strongest positive correlation (r=0.7, P<0.001), whereas the correlation between TET2 and NSUN5 was negative (r=−0.30, P<0.001). Weak correlations were observed between TRDMT1 and other regulators (YBX1, NSUN5, ALYREF, NSUN1, NSUN2, NSUN4, and ALKBH1) ( and Table S4). The PPI network analysis indicated that the 14 m5C regulators interacted with each other and TRDMT1 was one of the hub genes ().

Network analyses of m5C regulators

Univariate Cox regression analysis showed the prognostic values of 14 m5C regulators in patients with LUAD, and each regulator had a different prognostic value (). Among the m5C regulators, TRDMT1, NSUN1, NSUN4, NSUN7, and ALYREF were related to the prognosis of patients with LUAD (P<0.05). The interactions, connections, and prognostic values of the m5C regulators are depicted in the m5C regulatory network (). The strongest positive correlation was observed between TET2 and NSUN3 (r=0.70, P<0.001), while the strongest negative correlation was observed between TET2 and NSUN5 (r=−0.30, P<0.001) (Tables S5,S6).

Figure 2

Prognostic analysis of 14 m5C regulators and patterns of m5C methylation modification in TCGA-LUAD. (A) Prognostic analyses for 14 m5C regulators using a univariate Cox regression model. (B) The network of m5C regulators in LUAD. The lines linking regulators represent their interactions, and the thickness of the lines represents the correlation strength between regulators. (C) Kaplan-Meier curve analysis for patients with LUAD in clusters 1–3. (D,E) Functional annotation of DEGs among the 3 clusters using GO (D) and KEGG (E) analysis. (F) Heatmap and clinicopathological features of the three clusters classified by the m5C regulators. TCGA-LUAD, The Cancer Genome Atlas lung adenocarcinoma cohort; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Consensus clustering of m5C regulators identified 3 clusters with different clinical outcomes

To explore whether the expression levels of m5C regulators were associated with prognosis, consensus clustering analysis was applied to classify patients with LUAD in TCGA cohort into subgroups based on their consensus expression of m5C regulators. It was found that K=3 had optimal clustering stability to classify the patients with LUAD into 3 clusters, namely m5C clusters 1–3 (Figure S2A-S2H). Patients in m5C cluster 3 had a significantly poorer OS than patients in cluster 1 and cluster 2 (, P=0.032). Furthermore, to identify enriched functions and pathways among the clusters, GO and KEGG analyses were conducted based on the DEGs identified among the m5C clusters. The results indicated that the DEGs were enriched in various processes, including RNA transport, spliceosome, mitotic nuclear division, and chromosome segregation (). All significant (P<0.05) GO terms and KEGG pathways for the DEGs among 3 clusters are shown in Tables S7,S8. Furthermore, the associations among the 3 clusters, clinicopathological features, and expression levels of the 14 m5C regulators in the TCGA-LUAD cohort were evaluated. As shown in , the expression levels of the m5C regulators were higher in cluster 3, especially for ALYREF and YBX1.

Prognostic analysis of m5Csig in the TCGA-LUAD

As mentioned above, 5 m5C regulators, TRDMT1, NSUN1, NSUN4, NSUN7, and ALYREF, were associated with the prognosis of patients with LUAD according to the results of univariate Cox regression analysis (). The regulators TRDMT1, NSUN4, and NSUN7 act as protective factors, whereas NSUN1 and ALYREF are associated with risk of LUAD. The 5 m5C regulators were incorporated to build m5Csig according to the LASSO Cox regression algorithm. The regression coefficients of the 5 m5C regulatory factors are as follows: TRDMT1, −0.519056; NOP2, 0.166168; NSUN4, −0.376147; NSUN7, −0.246224; ALYREF, 0.163589. The patients with LUAD with complete clinical information (n=500) were classified into a high-risk group (n=250) and a low-risk group (n=250) according to the median m5Cscore, which was used as the cutoff point. As shown in , the expression levels of the risk-associated m5C regulators, NSUN1 and ALYREF, were upregulated in the high-risk group, and those of NSUN4, TRDMT1, and NSUN7 were downregulated in the high-risk group. With the increasing m5Cscore, the number of patients who died increased significantly (). The Kaplan-Meier curve revealed that the patients in the high-risk group had a significantly poorer OS than those in the low-risk group (P=6.578e-05, ). Principal component analysis (PCA) analysis showed that the high- and low-risk groups were stratified significantly in 2 different directions, indicating that the patients with LUAD in the high-risk group could be distinguished from those in the low-risk group (). The time-dependent ROC analysis indicated that the AUC values of m5Csig for 1-, 2-, and 3-year OS were 0.637, 0.615, and 0.658, respectively (). Next, univariate and multivariate Cox regression analyses were used to analyze whether the m5Cscore can be used as an independent prognostic factor. Univariate Cox regression analysis showed that T [hazard ratio (HR) =1.579, 95% confidence interval (CI): 1.296–1.923, P<0.001], N (HR =1.706, 95% CI: 1.405–2.072, P<0.001), M (HR =0.037, 95% CI: 1.038–3.272, P=0.037), stage (HR =1.577, 95% CI: 1.348–1.845, P<0.001) and m5Cscore (HR =2.800, 95% CI: 1.700–4.623, P<0.001) were significantly correlated with OS (). However, no significant correlation was found between age and gender and OS. Multivariate Cox regression analysis indicated that only m5Cscore (HR =2.263, 95% CI: 1.342–3.816, P=0.002) can be used as an independent prognostic factor for LUAD (). These results indicated that the m5Cscore has the potential to predict prognosis in LUAD patients.

Figure 3

Construction of m5Csig using 5 m5C regulators and the prognostic value of m5Csig in TCGA-LUAD. (A) A forest plot of univariate Cox regression identified 5 m5C regulators associated with overall survival. (B) The relationship between the expression profiles of the m5C regulators stratified by m5Cscore and the clinicopathological features of LUAD. (C,D) Survival status (C) and distribution (D) with increasing m5Cscore of patients with LUAD. (E) Kaplan-Meier survival curve of the high- and low-risk groups. (F) PCA of the TCGA cohort. (G) ROC curve and the AUC value of m5Csig for 1-, 2-, and 3-year overall survival. (H,I) Univariate (H) and multivariate (I) Cox regression analyses of clinicopathological features and m5Cscore. *P<0.05, **P<0.01. TCGA-LUAD, The Cancer Genome Atlas lung adenocarcinoma cohort; PCA, principal component analysis; ROC, receiver operating characteristic; AUC, area under the curve.

Validation of m5Csig in the GEO database

We validated m5Csig in the GSE72094 cohort. In total, 397 patients with complete clinical information were stratified into the high-risk group (n=198) and the low-risk group (n=199) using the median m5Cscore. As the m5Cscore increased, the number of deaths among the patients increased (). Patients in the low-risk group had a better OS than those in the high-risk group (P=1.58e-03, ). The PCA analysis suggested that patients were appropriately classified into high- and low-risk groups (). The ROC curves showed that the AUC values for 1-, 2-, and 3-year OS were 0.651, 0.615, and 0.59 respectively (). Univariate and multivariate Cox regression analyses also demonstrated that the m5Cscore can be used as an independent prognostic factor for LUAD patients ().

Figure 4

Validation of m5Csig in the GEO database. (A,B) Survival status (A) and distribution (B) of m5Cscore. (C) Kaplan-Meier survival analysis between the two groups. (D) PCA of the GEO cohort. (E) Time-dependent ROC and AUC based on the GEO data for 1-, 2-, and 3-year OS. (F-G) Univariate (F) and multivariate (G) Cox regression analyses of clinicopathological features and m5Cscore. GEO, Gene Expression Omnibus; PCA, principal component analysis; ROC, receiver operating characteristic; AUC, area under the curve; OS, overall survival.

Clinical characteristics between the high- and low-risk groups

A stratification analysis was performed to evaluate whether the m5Cscore could predict survival with the same clinical factor subgroup. Patients were stratified based on clinical parameters, such as age (£65/>65 years), gender (female/male), T (T1+2/T3+4), N (N0/N1–3), M (M0/M1), and stage (I+II/III+IV). The results showed that the m5Cscore could classify patients of the same stratum of age, gender, and early stage (T1–2, N0, M0 and stage I+II) into high- and low-risk groups (P<0.05). Patients in the high-risk group had a poorer OS than those in the low-risk group in each stratum (Figure S3A-S3L). We further analyzed the differences of clinical characteristics between the high- and low-m5Cscore groups and the difference of m5Cscore among different clinical characteristics. No significant distribution difference was found in terms of age (£65/>65 years) (P=0.15, Figure S4A,S4B), gender (female/male) (P=0.37, Figure S4C,S4D), stage I and stage II (P=0.19), stage I and stage III (P=0.33), and stage II and stage III (P=0.87) (Figure S4E,S4F). However, significant clinical differences were observed in terms of stage I and stage IV (P=0.043), stage II and stage IV (P=0.018), stage III and stage IV (P=0.023) (Figure S4E,S4F), ever smoking and never smoking (P=0.046, ), EGFR mutation group and EGFR wild group (P=4.2e-05, ), KRAS mutation group and KRAS wild group (P=0.0032, ), and TP53 mutation group and TP53 wild group (P=0.006, ).

Figure 5

Clinical characteristics and TMB between high- and low-risk groups. (A-D) The proportion and distribution of different clinical characteristics between high- and low-risk groups in GSE72094: ever smoking/never smoking (A,B), EGFR mutation/EGFR wild (C,D), KRAS mutation/KRAS wild (E,F), TP53 mutation/TP53 wild (G,H). (I) Difference in TMB between high- and low-risk groups. (J) Correlation scatter plots between TMB and m5Cscore. (K) Kaplan-Meier curves of OS for high- and low-TMB groups. (L) Kaplan-Meier curves of overall survival stratified by both TMB and m5Cscore. TMB, tumor mutation burden; OS, overall survival.

Tumor mutation burden in the high- and low- risk groups in the TCGA-LUAD database

The tumor mutation burden (TMB) quantification analyses indicated that the high-risk group correlated remarkably with a higher TMB (P<0.001, ). The m5Cscore and TMB also exhibited a significant positive correlation (R=0.24, P<0.001, ). There was no difference in OS between the high- and low-TMB groups (P=0.089, ). As shown in , when combined with the m5Cscore, there were significant survival differences among the 4 groups. The high-TMB/low-m5Cscore group had better survival than the high-TMB/high-m5Cscore group, and the low-TMB/high-m5Cscore group had the least favorable OS.

Expression of immune checkpoints and TME cell infiltration characteristics between the high- and low-risk groups in TCGA database

To determine the tumor immune infiltration characteristics, we evaluated the expression of 24 immune checkpoints between the high- and low-risk groups. We found a substantial difference in the expression of 24 immune checkpoints, among which LAG3, PDCD1 (PD-1), TNFRSF4, CD274 (PD-L1), CD276, TNFRSF8, TMIGD2, TNFRSF9, TNFSF4, TNFSF9, KIR3DL1, TNFRSF18, and CD70 were upregulated significantly in the high-risk group (). We further analyzed the proportion of immune cells between the high- and low- risk groups in the TCGA-LUAD database. Heterogeneity of LUAD was indicated by the different ratios of each cell type (). Furthermore, we compared the infiltration of immune cells between the high- and low- risk groups. As shown in , the high-risk group had a higher fraction of CD8 T cells, activated CD4 memory T cells, follicular helper T cells, resting NK cells, and M0 macrophages compared with those in the low-risk group (P<0.05). However, naive B cells, resting CD4 memory T cells, resting dendritic cells, and resting mast cells were markedly downregulated in the high-risk group (P<0.05).

Figure 6

Expression of immune checkpoints and TME cell infiltration characteristics between high- and low-risk groups in the TCGA database. (A) Expression of immune checkpoints between high- and low-risk groups. (B) Barplot of the distribution of 22 immune cells in TCGA-LUAD. (C) TME cell infiltration characteristics between high- and low-risk groups. *P<0.05, **P<0.01, ***P<0.001. TCGA-LUAD, The Cancer Genome Atlas lung adenocarcinoma cohort; TME, tumor microenvironment.

Construction of a prognostic nomogram for LUAD in the TCGA data

To establish a clinically applicable method to evaluate the prognosis of patients with LUAD, we constructed a prognostic nomogram by integrating clinical factors (age, gender, stage) with the m5Cscore (Figure S5A). Using the bootstrap method, calibration plots showed no significant deviation from the ideal for 1-, 3- and 5-year OS (Figure S5B). These results indicated that the prognostic nomogram could be used to predict the prognosis of patients with LUAD.

Discussion

Abnormalities of m5C modifications have been shown to influence RNA stability, gene expression, and protein synthesis, and thus have an essential role in various cellular, biological, and pathological processes (20-32). The RNA m5C modification and its regulators have been shown to be involved in the progression of various cancers, including hepatocellular carcinoma (33), bladder cancer (34), glioblastoma multiforme (35), breast cancer (36), and head and neck carcinoma (37), indicating that RNA m5C might play an important role in tumorigenesis and progression.

However, the biological functions and mechanism by which m5C modifications affect the TME were previously unknown.

In this study, we found that m5C regulators were significantly differently expressed in LUAD. By analyzing the expression profiles of m5C regulators, we identified 3 m5C clusters associated with different prognoses. Moreover, GO and KEGG function analysis indicated that DEGs among the 3 clusters were closely correlated with biological processes and signaling pathways, such as RNA transport, spliceosome, mitotic nuclear division, and chromosome segregation. We also identified 4 writers (TRDMT1, NSUN1, NSUN4, and NSUN7) and 1 eraser (ALYREF) were correlated with the prognosis of LUAD using LASSO Cox regression. Among the 5 m5C regulators, TRDMT1 mainly mediates tRNA stability and regulates cell metabolism of the m5C modification (30,38,39). Loss of TRDMT1 promoted homologous recombination and increased cellular sensitivity to DNA double-strand breaks (40). The NSUN1 protein (also known as NOP2) is a nucleolar-specific protein that plays a crucial role in RNA modification (41), cell cycle progression (41), chromatin organization (42), and HIV-1 latency (43). NSUN4, which forms a complex with MTERF4, is not essential in mitochondrial ribosome biogenesis and mitochondrial translation termination in conditional Nsun4 mouse knockout mutants (44,45). High expression of NSUN7 has been associated with shorter survival in low-grade gliomas (46), and a deletion mutation of NSUN7 has been associated with reduced sperm motility in asthenospermic men (47). The mRNA export adaptor, ALYREF, serves as a specific m5C-binding protein and functions in promoting mRNA export (48,49). An ALYREF-MYCN coactivator complex might be involved in neuroblastoma tumorigenesis (50).

An m5Csig was constructed, which divided patients with LUAD into high- and low-risk groups. Patients in the high-risk group had a significantly poorer OS than those in the low-risk group. Univariate and multivariate Cox regression analyses demonstrated that the m5Cscore was an independent prognostic factor for patients with LUAD. Accumulated evidence has demonstrated that patients overexpressing PD-1/PD-L1 and with a high TMB status are associated with an improved and durable ICB response (51-53). The TMB quantification analyses indicated that the high-risk group correlated markedly with a higher TMB. The m5Cscore and TMB also exhibited a significant positive correlation. The high-risk group displayed significantly higher expression levels of PD-1 and PD-L1 than the low-risk group. The above results demonstrated that LUAD with a high m5Cscore might show a better response to ICB therapy.

Accumulating evidence suggested that m5C is closely related to TME. Schoeler et al. demonstrated that TET enzymes control antibody production and shape the mutational landscape in germinal centre B cells. TET2 and TET3 guide the transition of germinal centre B cells to antibody-secreting plasma cells (54). Yue et al. revealed Tet2/3-deficiency in Treg cells leads to T cell activation and results in an activated phenotype and dysregulated expression of multiple Treg activation and phenotypic molecules in healthy mice (55). In our study, the CIBERSORT results showed that the high-risk group had stronger immune cell invasion compared with that of the low-risk group, for example, the numbers of CD8 T cells and activated CD4 memory T cells were significantly increased. These results suggested that the m5C regulators might be involved in the progression and prognosis of LUAD by modulating TIIC infiltration of the TME. Targeting m5C-related regulators might provide a novel way to improve the efficiency of ICB in LUAD.

However, there were several limitations to our study. First, this was a bioinformatic study based on a public database; therefore, further in vivo and in vitro experimental studies are needed to explore the potential effect and mechanism of m5C regulators in LUAD. Second, more potential m5C regulators have yet to be discovered. Last, the regulatory mechanism of m5C regulators in the TME was not determined, which requires further investigation to provide a better understanding.

Conclusions

In summary, we comprehensively analyzed the relationship between m5C methylation regulators and TME immune modulation. Based on the characteristics of m5C regulators, m5Csig was constructed to predict the prognosis of patients with LUAD, which might provide novel strategies for ICB therapy.

The article’s supplementary files as