Expression and prognostic significance of m6A-related genes in TP53-mutant non-small-cell lung cancer.

BACKGROUND: TP53 is an important tumor suppressor gene on human 17th chromosome with its mutations more than 60% in tumor cells. Lung cancer is the highest incidence malignancy in men around the world. N-6 methylase (m6A) is an enzyme that plays an important role in mRNA splicing, translation, and stabilization. However, its role in TP53-mutant non-small-cell lung cancer (NSCLC) remains unknown.
METHOD: First, we investigated 17 common m6A regulators' prognostic values in NSCLC. Then, after the establishment of risk signature, we explored the diagnostic value of m6A in TP53-mutant NSCLC. Finally, gene set enrichment analysis (GSEA), gene ontology (GO) enrichment analysis, and differential expression analysis were used to reveal the possible mechanism of m6A regulators affecting TP53-mutant NSCLC patients.
RESULTS: Study showed that nine m6A regulators (YTHDC2, METTL14, FTO, METTL16, YTHDF1, HNRNPA2B1, RBM15, KIAA1429, and WTAP) were expressed differently between TP53-mutant and wild-type NSCLC (p < 0.05); and ALKBH5 and HNRNPA2B1 were associated with the prognostic of TP53-mutant patients. After construction of the risk signature combined ALKBH5 and HNRNPA2B1, we divided patients with TP53 mutations into high- and low-risk groups, and there was a significant survival difference between two groups. Finally, 338 differentially expression genes (DEGs) were found between high- and low-risk groups. GO enrichment analysis, PPI network, and GSEA enrichment analysis showed that m6A may affect the immune environment in extracellular and change the stability of mRNA.
CONCLUSION: In conclusion, m6A regulators can be used as prognostic predictors in TP53-mutant patients.

INTRODUCTION

Lung cancer remains the deadliest malignancy in the world, and non‐small‐cell lung cancer accounts for 85% of all types of lung cancer individuals, which is characterized by high incidence and high mortality. Non‐small‐cell lung cancer can be mainly classified to three histological subtypes: lung squamous‐cell carcinoma (LUSC), lung adenocarcinoma (LUAD), and large‐cell lung cancer. Although great progress has been made in the detection and the treatment in NSCLC patients in recent years, the 5‐year survival remains for NSCLC remains only 16.6%.

TP53, also known as p53, is the most frequently mutated gene in NSCLC patients, even up to 80% in squamous‐cell carcinoma. It is one of the widely studied tumor suppressor genes, which names from the protein with the molecular weight of 53 kDa. P53 protein is an important regulator of cell growth, proliferation, and damage repair. DNA damage and oncogene activation can stimulate the acetylation and the activate of p53, to realize TP53 function as a cellular stress sensor. Moreover, TP53 is the most frequently mutated gene in human malignancies, and the tumor suppressive function of p53 protein in TP53‐mutant individuals was reversed compared with TP53‐wild type. In vivo experiments confirmed that Trp53 knockout mice had a higher risk of developing cancer.

People have discovered a reversible modification of RNA methylation that occurs at the sixth position of the RNA molecule adenine nitrogen atom (N‐6 methylation, m6A). M6A methylation is the most common post‐transcriptional modification in eukaryotes and plays an important role in mRNA metabolism and translation, as well as cell differentiation and embryonic development. It often found enriched in 3′‐UTR and near the termination codons of mRNA. M6A regulators can be classified into three types according to their functions, which are called “writers,” “erasers,” and “readers.” The function of the “writers” is modifying methyl to nucleotides, such enzymes include METTL3/14/16, RBM15/15B, WTAP, KIAA1429, and ZC3H13. METTL3, as a catalytic subunit, combines with METTL14 to form a hetero complex. And METTL3 was also found to play the role as a “reader” and located in the cytoplasm. WTAP can bind to the hetero complex and plays an important role in the recruitment of the hetero complex. METTL3/14‐WTAP complex can be induced into the nucleus by ZC3H13, thereby form the ZC3H13‐KIAA1429‐HAKAI complex in the nucleus to regulate the m6A process. RBM15/15B can also promote the methylation for certain RNAs. Methylation in mRNA caused by “writers” can be demethylated by “erasers” (including FTO and ALKBH5), which makes the process reversible. “Readers” is essentially a class of RNA‐binding proteins, which can be divided according to different domains, such as YTH domains (YTHDF1‐3 and YTHDC1/2), HNRNPs domains (HNRNPC and HNRNPB2A1), and some same RNA‐binding domains, respectively. , , These regulators can bind specifically to the methyl on mRNA. These genes build interaction network that worked by acting on m6A‐modified mRNAs.

So far, many studies have revealed the roles of m6A methylation modification in various tumors, particularly hepatocellular carcinoma, breast cancer, gastric cancer, and lung cancer. , , , However, the expression and prognostic significance in TP53‐mutant non‐small‐cell lung cancer are still unknown.

To further investigate the role of m6A modification in TP53‐mutant lung cancer patients, we conducted an in‐depth analysis of 17 m6A gene expression profiles in 469 lung cancer patients. In 469 patients, 233 individuals had a TP53‐mutant status. All the data were downloaded from the Cancer Genome Atlas (TCGA).

MATERIAL AND METHODS

Data set

All data were downloaded from TCGA database (https://portal.gdc.cancer.gov/), including 1026 NSCLC patients' RNA‐seq transcriptome profiling and 561 NSCLC patients' single nucleotide variation with corresponding clinical information.

Selection of m6A regulators

A total of 17 m6A regulators were finally identified in our analysis based on our search including eight “writers” (METTL3/14/16, WTAP, RBM15/15B, KIAA1429, and ZC3H13), seven “readers” (YTHDF1/2/3, YTHDC1/2, HNRNPA2B1, and HNRNPC), two “erasers” (FTO, ALKBH5). ,

Bioinformatics analysis

R 4.0.2 was applied to all analysis. First, we divided all samples into TP53‐mutant and TP53 wild‐type cohort. We compared expression levels between TP53‐mutant and TP53 wild‐type cohort. Univariate cox analysis was applied to independent prognostic analysis for 17 m6A‐related genes, hazard ratio (HR) value was used in determining protective or risk gene. Multivariate cox regression was used to construct an independent prognostic signature by regularizing and screening 17 genes in TP53‐mutant patients, and the risk score of screened genes was calculated. The calculation formula of risk score is . R package “Survival” and “glnmet” were applied to univariate and multivariate regression cox analyses in R version 4.0.2, respectively.

We then investigated the different pathways between high‐risk group and low‐risk group in TP53‐mutant cohort by gene set enrichment analysis (GSEA), and the GSEA test run 1000 times. Differentially expression genes (DEGs) were acquired between high‐ and low‐risk group in TP53‐mutant cohort. Ggplot2 package was applied to GO pathway enrichment analysis of these DEGs. A protein–protein interaction graph of these DEGs was constructed in STRING (http://string‐db.org/). MCODE plugin of Cytoscape software was applied to PPI visualization, and Bingo plugin was used to build GO pathways diagrams.

Statistics

Wilcoxon test was performed to compare m6A expression difference between mutant group and wild‐type group. Chi‐squared test was used to compare the difference of clinical features between different subgroup. Kaplan–Meier method was used for the analysis of overall survival. Univariate and multivariate cox regression analyses were applied to assess the relationship between prognostic and risk score. The missing data were deleted from the analysis. All the statistical analyses were performed in R version 4.0.2, and p value <0.05 was considered statistically significant.

RESULTS

Differential expression of m6A‐related genes in Lung cancer

Analyzing the expression of 17 m6A‐related genes in 233 TP53‐mutant samples and 236 TP53 wild‐type samples, we found that the expression of nine genes (YTHDC2, METTL14, FTO, METTL16, YTHDF1, HNRNPA2B1, RBM15, KIAA1429, and WTAP) were significantly different between two cohorts (Figure 1A,B). The correlation analysis showed that WTAP was significantly negatively correlated with YTHDC1, ZC3H13, METTL3, and METTL16 in the wild type, while the negatively correlated genes in mutant cohort were FTO and METTL3, YTHDC2, and HNRNPC. Most of the other genes were positively correlated in two cohorts (Figure 1C,D).

FIGURE 1

Differential expression of m6A regulators in TP53 wild‐type and mutant type and correlation in respective cohorts. (A). The heatmap of m6A regulators expression in every individual. (B). Violin plot of m6A regulators expression differential in wild‐type and mutant type. (C and D). Correlation in m6A regulators based on the Pearson correlation in wild and mutant cohorts. *p value <0.05, **p value <0.01, ***p value <0.005

Independent prognostic signatures building and comparison in wild‐type and mutant cohorts

We then constructed a prognostic model using univariate cox regression and multivariate cox analyses. Univariate cox regression was applied to screen prognostic‐associated genes. We found that HNRNPA2B1 gene has a good correlation with prognostics in patients without TP53 mutation. Results showed that ALKBH5 and HNRNPA2B1 lead to poor prognostic in TP53‐mutant patients. High expression of FTO and METTL14 meant different prognostic risk in patients between TP53‐mutant and wild‐type groups. Forest map showed that high expression of FTO and METTL14 was considered lower risk in the wild‐type patients, which was exactly the opposite in the mutant patients (Figure 2A,B).

FIGURE 2

(A and B) Univariate cox for m6A regulators in patients with and without TP53 mutant. (C) Mutant cases were divided into two groups by gene signature, and the Kaplan–Meier curves between two groups, with p value = 1.966e‐03. (D) ROC curve of gene signature to predict the prognostic in TP53‐mutant patients, AUC represents the sensitivity of prediction. (E) Multiple ROC curve of gene signature to compare with age, gender, stage, and TMN

Later, we constructed prognostic signatures of 17 genes in mutant cohorts using multivariate cox regression analysis to predict patients' prognostic risk score. Using this method, we established prognostic signatures containing ALKBH5 and HNRNPA2B1 in the mutant cohorts. Result of the multivariate cox regression and the coefficient value were shown in Table 1. Then, based on the risk score, patients were divided into the high‐ and low‐risk group. There was a significant difference between two groups in patients' overall survival (Figure 2C). The results of ROC curve also showed that it was feasible to evaluate the overall survival rate by risk signature in the TP53‐mutant patients (Figure 2D). Compared with clinicopathological features, risk score was also good predictors in prognosis (Figure 2E). To validate the signature in predicting patients' outcome, univariate and multivariate cox regression were applied to evaluate the accuracy of demography, clinicopathology, and risk score. In addition, the result showed that risk score performed well in predicting the prognosis of TP53‐mutant patients compared with clinicopathology features whether in univariate or multivariate cox regression analysis (Figure 3A,B). This result indicated that the independent prognostic signature constructed by ALKBH5 and HNRNPA2B1 had better predictive value than TNM in TP53‐mutant patients. Significant difference in gender was revealed by exploring the relationship between demographics and clinicopathology (Figure 3C). In the risk curves, the relationship of cases between the risk score and patients' survival were arranged in order of risk score from low to high, and we observed that alive patients became sparse as the risk score increases (Figure 3D,E). It can also be concluded that with the increment of HNRNPA2B1 and ALKBH5 expression, an increasing trend occurred in the risk score (Figure 3F). These findings confirmed the accuracy of the independent prognostic model.

TABLE 1

The coef value of the m6A regulators in risk signature.

N‐6 methylation regulators	Coef value	Hazard ratio
HNRNPA2B1	0.458124565	1.581105941
ALKBH5	0.514780837	1.673271743

FIGURE 3

(A and B) Verification of gene signatures by univariate cox regression and multivariate cox regression. (C) Associated between clinicopathological features and m6A regulators. (D and E) Risk score, survival time, and survival status of patients. (F) The relationship between expression of m6A regulators and prognostic risk. (G) Expression of HNRNPA2B1 and ALKBH5 in two groups. (H) Expression of all m6A regulators. HR, hazard ratio

Then, we explored the distinction of m6A regulators between high‐risk and low‐risk group. The results showed that the expression of ALKBH5 and HNRNPA2B1 in high‐risk group was higher than low‐risk group, and the difference was statistically significant (Figure 3G). Among them, more than a half of m6A regulators had a differential expression between two groups (Figure 3H).

M6A‐regulated signaling pathways and functional enrichment

Gene set enrichment analysis analysis was conducted to excavate the different signaling pathways and cell functional enrichment. As shown in Figure 4, biological functions related to cell proliferation and DNA synthesis were highly enriched in high‐risk group including cell cycle, spliceosome, folate‐involved one carbon metabolism, aminoacyl‐tRNA synthesis, RNA degradation, DNA replication, purine and pyrimidine metabolism, mismatch repair (MMR), and nucleotide excision repair (NER). IgA producing intestinal immune networks was observed to be silent in high‐risk group and enriched in low‐risk group (Figure 4). FDR <0.05 of signaling pathways and cell functions was used as inclusion criteria in GSEA analysis.

FIGURE 4

Performed GSEA analysis between two groups to seek different signaling pathways and cellular functions. It was revealed that cell functions in high‐risk patients is more associated with spliceosome, cell cycle, folate and one carbon pool metabolism, aminoacyl tRNA synthesis, RNA degradation, DNA replication, purine/pyrimidine metabolism, MMR and NER, etc. The IgA secreting intestinal immune network was found in low‐risk group, and that only this cellular function was enriched in low‐risk group and silenced in high‐risk group

Then, we searched for the differentially expression genes (DEGs) between high‐ and low‐risk groups. 338 DEGs were found by differential expression analysis. Heatmap and volcano plot were shown in Figure S1. The 338 DEGs were analyzed in BINGO plugin of cytoscape, and the functional enrichments of three modules of GO analysis (BP, CC, and MF) were shown. We found that most different pathways were located in extracellular region. Molecular functions dominated by these DEGs were more enriched in receptor binding, organic acid transmembrane transporter activity, IgE binding, symporter activity, etc. Biological processes were associated mainly with the regulation of response to external stimulus, cytolysis, humoral immune response, and response of fibrinolysis (Figure 5A). In addition, these molecular functions and biological process were also shown in the heatmap (Figure 5B).

FIGURE 5

(A) Tree modules (CC, BP, and MF) of GO analysis in cytoscape plugin BINGO. (B) Heatmap of GO functional enrichment analysis. (C) Protein–protein interaction network which accomplished in STRING. (D) Visualization of PPI networks in cytoscape. CC, cellular component; BP, biological process; MF, molecular functions; GO, Gene Ontology

We then put all the DEGs to STRING web tool (http://string‐db.org/) to analyze the protein–protein interaction network to further determine the molecular mechanism of the DEGS (Figure 5C). And the PPI network obtained from the String database is then imported in cytoscape for visualization. Green nodes represent the down‐regulated genes and red nodes represent the up‐regulated genes (Figure 5D).

Relationship between m6A regulators and prognostic in TP53‐mutant individuals

Finally, we assessed the relevance between all m6A regulators expression and overall survival in TP53‐mutant patients. All the cases were categorized into two groups by median expression value of respective gene, then Kaplan–Meier survival curves were plotted between high‐expression and low‐expression groups. According to the feedback information from the survival curves, we found four regulators including ALKBH5, METTL3, HNRNPA2B1, and YTHDC1 with their expression levels significantly correlated with survival (Figure 6).

FIGURE 6

The relationship between expression of each m6A regulator and overall survival in patients with TP53 mutations

Clustering and grouping prognoses for mutant and wild‐type cohorts

In order to explore a new m6A related classification of TP53‐mutant NSCLC, we performed the consensus clustering to all patients. The wild‐type and the mutant cohort were, respectively, divided into two clusters by consensus clustering, respectively, based on the expression profiles of all genes (Figure 7A,B). In the process of consensus clustering, we analyzed the possibility of clustering count (k value) from 2 to 9 (Figure 7C–H). Finally, k = 2 was applied to the mutant and the wild‐type cohort. The clinicopathological characteristics of the patients in each cluster were listed in Table 1.

FIGURE 7

(A–H), Cluster of respective cohorts

We then analyzed the clinicopathology and survival of two subgroups of each cohort. We found that two subgroups of wild type have a significant difference in T stage, while a more significant difference in M stage consists in mutant cohort. There was a significant difference in stage exists in both cohorts (Figure 8A,B). However, survival analysis, according to the results of both in the wild‐type cohort (p = 0.121) and mutant cohort (p = 0.089), survival situation between the subgroups were not significantly different (Figure 8C,D).

FIGURE 8

Survival and clinicopathology characteristics between two subgroups in TP53 wild‐type and mutant‐type correspondingly. (A and B), relationship between m6A regulator and clinicopathology; (C and D), Kaplan–Meier overall survival curves, elucidates the survival of the two subgroups. * differences are statistically significant

DISCUSSION

Lung cancer is the most common malignancy tumor among men and the second leading cause of death of malignancy in women that ranks only second to the breast cancer. In general, most patients were diagnosed with lung cancer with local or distant metastasis; hence, the average 5‐year survival rate for diagnosed lung cancer patients is as low as 15%. Therefore, systematic palliative is the most common clinical treatment for lung cancer, including radiotherapy, chemotherapy, and immunotherapy. As the most common mutated gene in human malignant tumors, TP53 gene plays a role in tumor formation, development, and treatment. About 50% patients have TP53 mutations, which are often missense mutations. Compared with wild‐type patients, TP53‐mutant individuals have poor prognostics and strong tolerance to chemoradiotherapy; hence, it can be used as a good prognostic biomarker for patients with NSCLC. P53 protein is the expression product of TP53 gene, which is often not detected in normal individuals. But the abnormal activation of p53 in cell can be caused while TP53 gene is in a mutated state, and the mutant protein can be endowed with functional acquisition. It is in this way that the TP53 gene participate in many stages with tumor development and mediates the treatment tolerance of tumors.

Chemical modification of RNA is considered to have important biological functions. So far, more than 160 RNA chemical modification have been found in eukaryotes of which N‐6 methylation that always as known as m6A is the most common. M6A was found functional in RNA transcript, RNA splice, and RNA degradation in mammals. Babieri I. showed that METTL3 and FTO could regulate the transcription of CEBP through its interaction with CEBP protein family and, thus, promotes the procession of AML. However, with the discovery of the first “eraser,” FTO, the m6A‐mediated methylation of the sixth nitrogen atom of RNA adenine has been recognized as a dynamic reversible regulation.

According to current studies, m6A regulators plays a role in variety of human tumors and can also be an important marker to predict prognostic in human malignancy. As the core gene of the “writers,” METTL3 plays a vital part in many tumors. In LUAD, METTL3 can augment EGFR expression and promote the cyclization of mRNA through eIF3. In addition, high expression of METTL3 can be considered as a signal of poor prognosis in lung adenocarcinoma, hepatocellular carcinoma, and gastric cancer. In lung squamous cell carcinoma, tumors with high expression of demethylase FTO show greater aggressivity and proliferation, as well greater apoptotic resistance (usually by acting on M2F1). Hence, FTO is clinically recognized as a prognostic factor for lung squamous cell carcinoma. ALKBH5 and HNRNPA2B1 were found to be predictors in TP53‐mutant patients in our study. ALKBH5 is a member of Alk protein family, with its dysregulation has been found in many tumors. ALKBH5 contains a DSBH domain, which can bind to the ATP domain of DDX3 gene, thus, affecting cell cycle, apoptosis, RNA degradation, and other cellular processes. In the present study, expression of ALKBH5 varied in different malignant tumors. ALKBH5 is often lower expression in colon cancer and pancreatic cancers, , while its higher expression is always a sign of poor prognosis in breast cancers and lung cancers. , We focused on expounding the role of ALKBH5 in NSCLC. Evidence suggested that reduced levels of n‐6 methylation of FOXM1 can inhibit proliferation and invasion phenotypes of lung cancer cells. This reduction occurred in ALKBH5 knock‐down individuals. HNRNPA2B1 is closely related to a function called m6A‐swicth. It has a HNRNPs domain which has specific structures, thus, to regulate RNA–protein interactions. In addition, it can regulate the splicing of transcriptional exons and the formation of pre‐miRNA, which is similar to METTL3—m6A “readers”.

In our study, we were devoting to construct a signature that contains m6A regulators that were expected to predict the prognostic in TP53‐mutant non‐small‐cell lung cancer patients. We first analyzed the differences of m6A regulators' expression in TP53‐mutant and wild‐type NSCLC patients and found that more than a half of m6A regulators expression were different between two groups. The following univariate cox regression analysis revealed that m6A regulators predicted different prognostic risk with different status of TP53 gene. For example, high expression of ALKBH5 in TP53‐mutant patients meant a poor outcome, but with no significant effect in TP53‐mutant individuals. We then established the prognostic signature in patients with TP53 mutant. Multivariate cox regression showed that a signature consisting ALKBH5 and HNRNPA2B1 could categorize the patients into high‐ and low‐risk groups. Kaplan–Meier survival curve indicates that there are significant differences between two groups obtained by the signature. The validity of the prediction was verified by univariate independent prognostic and multivariate independent prognostic analysis. The forest map showed that risk score was a better predictor of patients' prognostic than stage and TNM. The risk curve also confirmed the conclusion. The majority of m6A regulators in high‐risk group has a significantly higher expression. GSEA enrichment analysis was carried out between high‐ and low‐risk groups next to analysis possible mechanism. It was found that the enriched pathways in high‐risk group were more related to cell proliferation and gene expression. The conclusion predicts the possible impact of m6A regulators in patients with TP53 mutations. And the DEGs between groups were identified for further functional exploration. These DEGs were put in GO enrichment analysis and visualized. The results showed that these DEGs were significantly enriched in functions as receptor binding, cytolysis, and so on, which mostly occurred in extracellular region and plasma membrane, which meant m6A may affect the microenvironment of TP53‐mutant cancer cells. Finally, we used the method of the consensus clustering to classify the TP53‐mutant and wild‐type cohorts into any subgroups, respectively, in an attempt to find a better prediction method for tumor prognosis.

In our study, we systematically analyzed the potential value and possible mechanism of m6A regulators in TP53‐mutant NSCLC patients. However, there were still some limitations in this study. First, all the data were downloaded from the TCGA database and should be validated in lung cancer patients; second, all the cases were from the United States, which could lead to regional and racial biases in the results; third, the exact mechanism remains unclear.

In short, m6A regulators may provide new insights and explore new target for the diagnosis and treatment of patients with TP53 mutations in NSCLC patients.

CONFLICT OF INTEREST

The authors declare that there is no relative conflict of interests.

AUTHOR CONTRIBUTIONS

Zhuochen Zhao, Junhu Wan, and Manman Guo were jointly responsible for the data acquisition and analysis; Yangxia Wang, Zhengwu Yang, and Zhuofang Li involved in the production of the Figures and Tables; Zhuochen Zhao and Junhu Wan conceived and wrote the article; Fuyou Zhou involved in the revision of the article. Liang Ming approved the final article.

CONSENT FOR PUBLICATION

Written informed consent for publication was obtained from all participants.

Fig S1

Click here for additional data file.

Tab S1

Click here for additional data file.