Identification of a novel glycolysis-related prognostic signature for predicting prognosis and tumor microenvironment of lung adenocarcinoma

doi:10.21203/rs.3.rs-2375858/v1

Download PDF

Research Article

Identification of a novel glycolysis-related prognostic signature for predicting prognosis and tumor microenvironment of lung adenocarcinoma

https://doi.org/10.21203/rs.3.rs-2375858/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Glycolysis affects tumor growth, invasion, treatment resistance, and the tumor microenvironment. In this study, we aimed to construct a glycolysis-related prognostic signature for lung adenocarcinoma (LUAD) and analyze its relationship with the tumor microenvironment (TME).

Methods: We analyzed the data of a training set from The Cancer Genome Atlas (TCGA) database and four validation cohorts from the Gene Expression Omnibus (GEO) databases which included 1,689 patients with LUAD. The genetic and transcriptional alterations of glycolysis-related genes (GRGs) were investigated, and evaluated their prognostic value in LUAD patients. The glycolysis-related patterns were identified using consensus unsupervised clustering analysis. A glycolysis-related prognostic signature was then established using the least absolute shrinkage and selection operator and Cox regression analysis. Gene set variation, clinical relevance, and TME analyses were conducted to explain the biological functions of glycolysis regulators and their performance in prognostic prediction.

Results: The glycolysis-related patterns were identified based on 18 GRGs, with significant differences in survival between the patterns. By screening of differentially expressed genes between patterns, a novel glycolysis-related prognostic signature was final obtained, including ALDOA, FOSL2, PDE6D, PPARD, and RASAL2, with robust and accurate prognostic performance. The high-, and low risk groups, distinguished by the signature formula, presented a significant difference in overall survival, TME, and immunotherapy response in LUAD patients. The potential roles of the glycolysis-related prognostic signature were confirmed using the validation dataset.

Conclusions: This study identified five glycolysis-related prognostic genes that effectively predicted the survival of LUAD. To a certain extent, the newly identified gene signature was related to the TME, especially immune cell infiltration. These findings provide potential biomarkers and therapeutic targets for LUAD.

glycolysis

lung adenocarcinoma

tumor microenvironment

immunotherapy efficacy

prognostic prediction

Lung cancer is the leading cause of cancer death worldwide [1], non-small cell lung cancer (NSCLC) occupies 85% of lung cancer cases, and lung adenocarcinoma (LUAD) is one of the major pathological types of NSCLC [2]. Although breakthroughs have been made in the molecular targeted therapy and immunotherapy of LUAD, the long-term survival rate is still non-ideal [3, 4]. Therefore, an urgent need exists to investigate the pathogenic mechanisms underlying lung cancer and establish a reliable predictive model for the prognosis of patients diagnosed with LUAD.

The deregulation of cellular energetics or metabolic reprogramming is a key characteristic of tumor cells, closely related to tumor occurrence, progression, and drug resistance [5]. Normal cells rely on glycolysis, instead of oxidative phosphorylation, to provide energy during hypoxia. However, in tumor cells, glycolysis tends to occur even in aerobic environments. This phenomenon of tumors is known as the “Warburg effect”, which has been found to play significant role in tumor invasion, immune escape, angiogenesis, autophagy, and tumor microenvironment (TME) [6–9]. For example, the lactic acid produced by glycolysis accumulates extracellular, which reduces tissue pH and affects the TME in turn promotes to tumor progression and treatment-resistant [10–12]. In addition, through concentration gradient-dependent inhibition, excessive lactic acid in the TME inhibits T cells from releasing lactic acid, thereby affecting the proliferation and function of T cells [13].

Little is known regarding glycolysis-related genes (GRGs) with a predictive value and the relationship between glycolysis and TME in LUAD. In this study, we aimed to clarify the potential role of GRGs in LUAD. By comprehensively analyzing copy number variations (CNVs), mutation annotations, methylation levels, and gene expression profiles from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases, we identified potential pathways related to glycolysis of LUAD. Subsequently, the robust and effective prognostic signature were constructed on basis of differentially expressed genes (DEGs) identified based on the two glycolysis-related patterns, and evaluated its underlying connection with the tumor microenvironment (TME) and response to immunotherapy. Our results provide new insights for further research on glycolysis of tumor and can inform personalized LUAD management.

2.1. Data acquisition and preprocessing

The transcriptome profiling data, gene methylation data and corresponding clinical data of LUAD were acquired from the TCGA database and GEO database, 62 GRGs were retrieved from MSigDB database for subsequent bioinformatics analysis. The TCGA-LUAD dataset including 535 tumor and 59 tumor-adjacent samples was selected as the training cohort. For the GEO datasets, the GSE13213 [14] and GSE68465 [15] cohorts were used for internal validation, whereas the GSE31210 [16] and GSE72094 [17] cohorts for further external validation. Supplementary Table S1 list the detailed clinical characteristics of the whole 1,689 patients with LUAD in this study. Normalized expression values were log₂-transformed and scaled before use in model validations. Using the “combat” algorithm in the “SVA” package of R software, we corrected the batch effect between TCGA and GEO datasets. For analyses involving clinical data, samples associated with a survival time less than 30 days or unknown were removed.

2.2. Identification of differential expression and genetic alterations in GRGs

The “limma” R package was used to screen DEGs with a threshold p-value < 0.05 between tumor and normal lung tissues from TCGA-LUAD dataset. Univariate Cox analysis was further used to select GRGs with prognostic value from the DEGs. We determined CNVs using GISTIC 2.0, including amplification (segment mean > 0.2), diploid (–0.2 < segment mean < 0.2), and deletion (segment mean <–0.2). Each type of CNVs for each glycolysis-related prognostic DEG was visualized with the R package “RCircos.” The “maftools” R package was used to determine the genetic mutation frequency and types of glycolysis-related prognostic DEGs in TCGA-LUAD samples.

2.3. Identification of GRGs associated with survival and hierarchical clustering

Unsupervised consistent clustering analysis was performed to classify the LUAD samples from the training and internal validation cohorts into different patterns based on the expression levels of the glycolysis-related prognostic gene set. The number of clusters was determined using the R package “ConensusClusterPlus” with 100 replicates and 80% item resampling to verify the stability of the clusters. Principal component analysis (PCA) was used to determine the ability of molecular clusters to discriminate between patients in subgroups as a dimensionality reduction method, and the results were visualized using the “FactoMineR” algorithm.

2.4. Construction and validation of a prognostic model

To construct a glycolysis-related risk model with better prognostic value, we focused on the DEGs identified between two glycolysis-related patterns. The DEGs were filtered with a threshold adjusted p-value < 0.001 and │log₂(fold change)│ > 1. The least absolute shrinkage and selection operator (LASSO) regression method with 10-fold cross-validation was used to narrow down the prognosis-related DEGs applying the R package “glmnet.” Multivariate Cox regression analysis was performed to establish a signature for evaluating the relationships between the DEGs and the survival of LUAD patients. Finally, the risk score based on the glycolysis-related risk model was constructed with the following formula: Risk score = Ʃi Expression (mRNA_i) × Coefficient (mRNA_i).

2.5. Evaluation of the risk model of GRGs as an independent prognostic factor in LUAD

Both the training and validation cohorts were classified into high- and low-risk groups, with the median risk score used as the cut-off point. The predictive value of the prognostic model was assessed using a time-related receiver operating characteristic (ROC) curve. Heatmaps were used to compare and visualize the distribution of clinicopathological characteristics between the risk cohorts. Univariate and multivariate Cox regression analysis were performed to test the prognostic independence of the model. The risk scores of each patient were further combined with clinical characteristics to construct a nomogram through the “rms” R package.

2.6. Enrichment analysis

Gene Ontology (GO) was used to annotate the biological processes, molecular functions, and cellular components of genes. Gene pathways were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG). Gene set variation analysis (GSVA) was performed to assess differences in biological pathways between patterns using the R packages “GSEABase”, “GSVA”, and “pheatmap”. The annotated gene set “c2.cp.kegg.v7.0. symbols.gmt” was used as the reference gene set, with FDR < 0.05 as the screening threshold.

2.7. Assessing the response to immunotherapy and drug sensitivity analysis

At present, it is widely recognized that immune checkpoints are the major mechanism of immune resistance, especially against T cells specific for tumor antigens [18]. To analyze the response to immunotherapy in distinct risk cohorts, we assessed the potential response to immune checkpoint inhibitor (ICI) therapy according to Tumor Immune Dysfunction and Exclusion (TIDE) and the immunophenoscore (IPS) [19]. The Estimation of STromal and Immune cells in MAlignant Tumor tissues using expression data (ESTIMATE) algorithm was applied to determine the components of stromal and immune cells according to the gene expression characteristics of the tumor samples. Many studies have reported a potential link between RNA modification and immune cells [20]. Hence, we used the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm to evaluate the abundance of diverse immune cells and immune functions in the samples of different clusters. Moreover, the R package “oncoPredict” was used to analyze the expression profile of TCGA-LUAD cohort in drug sensitivity analysis [21].

2.8. Statistical analysis

All statistical analyses were performed using R v4.2.2 software. Spearman correlation analysis was employed to calculate the correlation coefficient between variables. The median risk value was used as a cut-off value to divide the patients into high-risk and low-risk groups. The Kaplan–Meier (K–M) method was used to compare the overall survival (OS) between the groups. p-values < 0.05 (*), 0.01 (**), and 0.001 (***) were considered statistically significant.

3.1. Genetic variation of glycolysis-related prognostic genes in LUAD

In this study, we evaluated 62 GRGs in LUAD, and 28 DEGs were screened between LUAD and normal lung tissues from the TCGA-LUAD dataset (Supplementary Figure S1A, B). Univariate Cox survival analysis further identified 18 glycolysis-related prognostic DEGs from the above genes (Supplementary Table S2). As shown in Fig. 1A, we found a higher frequency of CNV amplification in GPI, ALDH3B2, LDHB, and ALDOA, whereas extensive CNV deletions were identified in DLAT, PGM1, ALDH2, and ENO3. Figure 1B shows the locations of LUAD-associated CNVs in glycolysis-related prognostic DEGs on chromosomes. As shown in Supplementary Figure S1C, CNV amplification was associated with increased ALDH1B1, DLAT, LDHB, and LDHA, et al. expression levels, which verified that CNVs play a potential role in regulating the expression of glycolysis-related prognostic DEGs.

We then analyzed the correlation between the methylation and mRNA expression of glycolysis-related prognostic DEGs (Supplementary Figure S2), and found that gene methylation was associated with decreased ALDH2, ENO3, FBP1, and PGM1, et al. expression levels. Additionally, mutations in glycolysis-related prognostic DEGs were identified in only 71 of the 539 samples (13.2%) in TCGA-LUAD cohort (Fig. 1C). Considering that the frequency of CNVs and somatic mutations in glycolysis-related prognostic genes were low overall, this indicated that GRGs may be conserved in LUAD. Comparison of glycolysis-related prognostic genes expression in LUAD and normal lung tissue samples demonstrated that the majority of genes were upregulated in tumor samples (Fig. 1D).

The protein expression of GRGs was further validated using the immunohistochemistry results from the Human Protein Atlas platform. As shown in Fig. 1E and Supplementary Figure S3, the protein expression levels of ALDOA, ENO1, GAPDH, LDHA, and PFKP were higher in tumor tissues than in normal tissues. The PCA plot showed that the expression profiles of the 18 GRGs could clearly distinguish between normal and LUAD samples (Fig. 1F). Collectively, these results indicated that the different expression profiles and genetic variation landscapes of GRGs could play an important role in LUAD tumorigenesis and progression.

3.2. Evaluation of glycolysis-related patterns based on 18 regulators

K–M survival analysis from combined samples of the TCGA, GSE68465, and GSE13213 cohorts suggested that LUAD patients with higher expression of ADH1B, ALDH2, ENO3, and FBP1 showed improved clinical outcomes, whereas higher expression of the rest genes was associated with worse OS (Supplementary Figure S4). The gene network of GRGs is shown in Fig. 2A, depicting their interactions, correlation of mRNA expression, and prognostic value for LUAD patients. These data revealed that most GRGs can be considered as risk factors for OS in TCGA-LUAD cohort.

Consensus cluster analysis was performed to classify the patients based on the expression of 18 GRGs that were significantly associated with OS. Two glycolysis-related patterns were identified, with 666 and 393 patients in patterns A and B, respectively (Fig. 2B–C). PCA suggested an obvious intergroup distribution between the two patterns (Fig. 2D). The box plot in Supplementary Figure S5A shows the expression levels of the 18 glycolysis-related prognostic genes in the two patterns, demonstrating that pattern A presented higher expression of GRGs. And there was significant difference in OS between the two patterns, pattern B showed a trend toward superior survival (Fig. 2E). GSVA showed that compared to glycolysis-related pattern B, pattern A was significantly enriched in certain KEGG pathways associated with metabolism and base repair, such as pentose-phosphate pathway, pyrimidine and purine metabolism, mismatch repair, and nucleotide excision repair (Fig. 2F and Supplementary Table S3).

3.3. TME cell infiltration characteristics of glycolysis-related patterns

Based on ssGSEA, pattern B had a higher abundance of most types of immune cells than pattern A (Fig. 3A). Similar results were obtained in assessing the immune functions of the two clusters (Supplementary Figure S5B), which may partly explain the better OS in cluster B than in cluster A. Evaluation based on the ESTIMATE algorithm revealed that the stromal score, immune score, and ESTIMATE score were all decreased in the pattern A (Fig. 3B). Subsequently, we evaluated the response to ICI therapy using the TIDE algorithm. As speculated, immunotherapy was more likely to be effective in patients in pattern B than in those in pattern A according to the TIDE score (Fig. 3C). Spearman correlation analysis and the TIMER algorithm [22] showed that the expression levels of the GRGs ADH1B, ENO3, and PGAM1 significantly correlated with the infiltration of many types of immune cells (Fig. 3D and Supplementary Figure S6).

3.4. Construction of the glycolysis-related clusters

By considering the 5582 DEGs between two glycolysis-related patterns from the combined datasets, we conducted GO enrichment, and KEGG pathway analyses to explore the potential biological functions of these genes. As shown in Fig. 4A and Supplementary Figure S7A, “regulation of mitotic cell cycle” was the most enriched terms among the biological process categories, and “actin cytoskeleton” and “DNA-binding transcription factor binding” were the most enriched terms among the cellular component and molecular function categories, respectively. Pathways of neurodegeneration was identified to be the most enriched among the KEGG pathways of the DEGs (Fig. 4B and Supplementary Figure S7B). Univariate Cox analysis was used to identify 471 prognostic related genes from the above DEGs (Supplementary Table S4) and consensus cluster analysis was performed to identify GRG clusters (Supplementary Figure S7C, D). Comparison of 18 glycolysis-related prognostic genes expression in glycolysis-related clusters demonstrated that the majority of genes were significantly upregulated in cluster A (Supplementary Figure S7E). The K–M curve showed significant alterations in survival outcomes between the two glycolysis-related clusters, the cluster B presented a clear OS advantage, whereas the cluster A had a higher risk of mortality (Supplementary Figure S7F).

3.5. Construction of the glycolysis-related prognostic signature

To promote the clinical application of glycolysis-related patterns, LASSO regression was adapted to filter 471 candidate genes down to five genes in combined datasets based on the minimum criteria used to construct the glycolysis-related risk model (Fig. 4C, D). The formula was as follows: risk score = (expression level of ALDOA × 0.173) + (expression level of FOSL2 × 0.211) + (expression level of PDE6D × 0.175) + (expression level of PPARD × 0.235) + (expression level of RASAL2 × 0.160). The risk score of each sample in the TCGA, GSE68465, and GSE13213 cohorts was calculated and patients were divided into two groups according to the median risk score (Supplementary Table S5).

As shown by the K–M curves, patients in the high-risk group had significantly worse OS than those in the low-risk group (p < 0.01, Figs. 5A–C); analogously, as the risk score increased, more patients died (Supplementary Figure S8A–C). In the training cohort, the AUC values of the present risk model were 0.646, 0.636, and 0.592 for the 1-, 3-, and 5- year prognoses, respectively (Fig. 5D), with similar results observed in the internal and external validation groups (Figs. 5E, F). Univariate and multivariate Cox analyses indicated that the glycolysis-related signature was an independent prognostic factor for OS in TCGA-LUAD cohort (Fig. 5G), which consistent with the results observed in the internal and external testing dataset (Fig. 5I–L).

We further investigated the association among the glycolysis-related patterns, gene clusters, risk score, and survival state, which is visualized as a Sankey diagram in Supplementary Figure S8D. As shown in Supplementary Figure S8E, F, pattern A and cluster A exhibited a significantly higher risk score than pattern B and cluster B, respectively. Considering that pattern A and cluster A were associated with a worse survival trend, this result again demonstrated the prognostic value of the glycolysis-related risk model.

The distribution of clinicopathological features between the high- and low-risk groups from combination of training and internal validation cohorts are presented as a heatmap in Fig. 6A and as a bar chart in Fig. 6B–D. These results suggest that patients with a higher stage LUAD were significantly associated with a higher risk score. K–M curve showed that the glycolysis-related risk model exhibited prognostic predictive power in various subgroups with different clinical characteristics (Fig. 6E–G), with similar results observed in the external validation group (Supplementary Figure S9). To facilitate the utilization of our risk model, we constructed a nomogram using risk scores combined with clinical characteristics, as shown in Supplementary Figure S10A, which further verified the reliability and accuracy of this nomogram in predicting 1-, 3-, and 5-year OS (Supplementary Figure S10B). These results demonstrated that the glycolysis-related risk model is stable and can precisely predict the clinical outcomes of LUAD patients.

To uncover the potential biological characteristics between high-, and low-risk groups, we performed GSVA enrichment analysis. In addition to pathways related to base repair, the high-risk group was also significantly enriched in various tumor pathways such as small cell lung cancer and renal cell carcinoma, compared with the low-risk group (Fig. 7A and Supplementary Table S6). The activation of these pathways could be a potential mechanism underlying the poor survival in high-risk groups.

3.6. Correlation of the glycolysis-related prognostic signature with the TME and immunotherapy response

The difference between the two glycolysis-related risk groups in carcinogenesis pathways prompted further analysis of TME cell infiltration. As shown in Fig. 7B and Supplementary Figure S10C, the glycolysis-related risk model was associated with the abundance of most immune cells and immune function. In addition, we observed that patients in the low-risk group had higher immune score (Supplementary Figure S10D). We further analyzed the underlying links between the risk score and tumor mutation burden (TMB) in the TCGA-LUAD cohort. the risk score was positively correlated with the TMB (R = 0.3, p < 0.001) (Supplementary Figure S10E, F). The TMB landscape suggested that the most common mutations were found in TP53, TNN, and MUC16 (Supplementary Figure S10G, H).

Additionally, patients in the low-risk group had a lower TIDE score (Fig. 7C) and higher IPS (Fig. 7D-G) were observed, which suggested that LUAD patients in the low-risk group may benefit more from immunotherapy. To determine the effectiveness of the risk score for predicting the drug treatment response in patients with LUAD, we evaluated the response to chemotherapeutic and small-molecule drugs between the two groups. In the high-risk group, 67 drugs with markedly lower half-maximal inhibitory concentration (IC50) values were observed; concomitantly, 19 drugs were observed in the low-risk group (Supplementary Table S7).

Metabolic reprogramming has been shown to play important roles in tumor cell growth, proliferation, and survival [5]. Recent research has focused on glycolytic processes in cell invasion, autophagy, and treatment resistance [6–9]. Despite extensive research on prognostic signatures in LUAD, studies involving GRG with a predictive significance are insufficient.

In the present study, we initially revealed the transcriptional expression patterns of GRGs, and the prognosis was further analyzed. Two distinct molecular patterns were identified based on 18 prognostic GRGs. By comparing with patients with pattern A, patients with pattern B had better OS. According to GSVA enrichment analysis, pattern A was significantly enriched in base repair and metabolism related pathways. The level of immune cell infiltration and immune function also different greatly between the two patterns. Subsequently, a novel glycolysis-related signature based on 5 prognostic related DEGs, including ALDOA, FOSL2, PDE6D, PPARD, and RASAL2, was generated using Lasso Cox hazards regression analysis in the testing cohort and internal validation cohort, then calculated the risk score based on the signature formula.

Most of these genes have been studied in a variety of tumors. GRG ALDOA encodes glycolytic enzyme aldolase A. It is reported that the established effect of ALDOA including epithelial–mesenchymal transition [23] and multiple growth factor signaling pathways crosslink with alternative metabolism, such as the EGFR/MAPK and PI3K/AKT-mTOR pathways [24, 25]. Moreover, ALDOA phosphorylation was validated to promote glycolysis and cell proliferation in the hepatocellular carcinoma [26]. Recent studies have shown that ALDOA is a key enzyme involved in lung cancer metabolic reprogramming and metastasis. Overexpression of ALDOA is related to various clinical parameters and increase migration and invasion of lung cancer cell lines in vitro and formation of metastatic lung cancer foci in vivo.

Fos-like antigen 2 (FOSL2/FRA-2) is a member of the activator protein 1 transcription factor family. As a transcriptional switch in the lung tumor microenvironment, the activation of Wnt/β-catenin pathway-mediated FOSL2 drive transformation of gene regulatory from M1-like tumor-associated macrophages (TAMs) to M2-like TAMs, thereby promoting lung cancer progression and metastasis [27]. Furthermore, Wan et al., reported that the enhancement of FOSL2 in breast cancer-associated fibroblasts was significantly associated with clinical progression by promoting angiogenesis [28].

PDE6D encodes phosphodiesterase 6δ, a RAS chaperone protein, and was found to mediate the antegrade trafficking of prenylated KRAS to the cell membrane, where KRAS is activated [29, 30]. Studies suggested that targeting PDE6D-KRAS interaction is an effective approach to tackle both KRAS wild type and mutated cancer types [31–33]. Moreover, PDE6D has been found to be significantly expressed in cancer cells and tissues and was considered as potential novel biomarker for colorectal cancer [34]. However, there is no evidence yet to suggest an underlying role of the trafficking chaperone PDE6D in lung cancer.

Peroxisome proliferator-activated receptor delta (PPARD) is a ligand-dependent nuclear receptor that functions as a transcription factor has pleiotropic effects on cell homeostasis [35]. The expression of PPARD is upregulated in many cancers, such as breast [36], gastric[37], and colon [38]cancer and associated with significantly reduced metastasis-free survival. Zuo et al. [39] reported that PPARD expression in cancer cells drastically affected epithelial-mesenchymal transition, migration, and invasion, further underscoring its necessity for metastasis. These founds demonstrate that PPARD is an important molecular target in pan-cancer.

RAS protein activator like 2 (RASAL2) is a member of the family of RAS GTPase-activating proteins. The RASAL2 protein negatively regulates the RAS signaling pathway by catalyzing the hydrolysis of RAS-GTP to RAS-GDP in many cellular activities and acts as a vital regulator of the RAS signaling pathway. However, its tumor-suppressive or oncogenic roles in cancer development remain controversial. Studies indicate that the biological function of RASAL2 is influenced by cellular context to influence its pro or anti-oncogenic activity in human cancers and different signaling pathways of RAS signaling pathway, PI3K/AKT/mTOR signaling pathway, and NF-κB pathway may account for the different biological outcomes of RASAL2 activity [40–44].

The GSVA revealed that the high-risk group was mainly enriched in proliferation and survival pathways associated with cell cycle, DNA replication, and base repair. Patients in the low- and high- risk score groups showed significant differences in TME, TMB, and drug sensitivity. Univariate and multivariate Cox analyses confirmed that the novel glycolysis-related signature was an independent prognostic indicator for LUAD. ROC curves further showed the efficacy of the glycolysis-related risk score in predicting survival rate. Similar results were obtained from the independent validation cohort, demonstrating that the glycolysis-related signature may be a robust biomarker for evaluating the prognosis of LUAD. Finally, we established a nomogram for clinical use by combining the glycolysis-related signature with several clinicopathological features.

Despite recent advances in immunotherapy for LUAD, choosing the appropriate therapeutic regimen for patients remains a clinical challenge [45]. In addition, some patients do not obtain effective benefits from immunotherapy [46–48] and some patients experience severe side effects during therapy [49]. Therefore, it is critical to explore a novel method to guide the individualized and precise treatment of patients with LUAD. We speculated that patients stratified using the glycolysis-related signature would show distinct immunotherapeutic responses, which was investigated based on the associations among glycolysis-related score groups, TME factors, and the TIDE score.

The highly acidic microenvironment produced by tumor glycolysis may affect the infiltration of immune cells to varying degrees, eventually leading to tumor progression and immune escape [50]. In this study, we found that the infiltration of CD4⁺ T cells, CD8⁺ T cells, DCs, and neutrophils differed significantly between high- and low-risk groups, demonstrated that glycolysis affects the composition of immune cell populations in TME of LUAD. It was reported that TIDE could be an accurate biomarker to predict the immunotherapeutic effects in non-small cell lung cancer, which is negatively associated with the efficacy of ICIs [19]. Another recent study demonstrated the clinical application of TIDE in predicting and evaluating the immunotherapeutic response [51]. In our study, compared with patients in the high-risk group, those with lower risk scores had lower TIDE scores, suggesting that patients in the low-risk group may obtain more clinical benefits from immunotherapy, which was verified in the IPS analysis. Collectively, these findings highlight a crucial role of glycolysis in mediating the clinical response to ICIs therapy by affecting immune cell infiltration. Thus, the glycolysis-based signature is a potential biomarker for assessing immunotherapeutic response and tailoring individualized treatment for patients with LUAD.

This study also has some limitations that should be acknowledged. First, despite the large scale of sample size, the data included in this study were derived from public databases, resulting in selection bias, thus affecting the generalizability of the results. Second, the potential functional mechanisms of the glycolysis-based risk score have not been fully verified, which requires follow-up functional experiments at the molecular level both in vivo and in vitro. Finally, there are missing data for clinical variables such as targeted therapy and treatment after recurrence, which require improvement, and more clinical variables should be integrated to further validate the clinical value of the glycolysis-related prognostic signature.

In conclusion, the present study systemically analyzed the genetic characteristics and prognostic value of GRGs in LUAD, and emphasized the vital role of the glycolysis in shaping the heterogeneity and complexity of the TME. We constructed a novel prognostic model based on glycolysis-related signature to predict survival and guide treatment decisions. The findings of this study provide insights into the glycolysis-regulated TME in LUAD and identify potential immunotherapeutic targets as a new approach for the clinical individualized treatment of LUAD.

Funding information

Not Applicable.

Acknowledgements

We thank the builders of the Harmonized Cancer Datasets, the NCI and the National Human Genome Research Institute in the US; as well as the GSE96058 dataset creators Tomida S, Takeuchi T, and Shimada Y, et al., from Aichi Cancer Center, Japan; GSE68465 dataset creators Shedden K, Taylor JM, and Enkemann SA, et al., from NCI, USA; GSE31210 dataset creators Hirokazu O, Takashi K, and Yuko I, et al., from National Cancer Center Research Institute, Japan; GSE72094 dataset creators Schabath M, Welsh E, and Fulp W, et al., from Moffitt Cancer Center, USA. We are also grateful for the contributions of all patients and their families who participated in this study.

Disclosure

All authors have declared no conflicted of interest.

Ethics approval and consent to participate

This study was approved by the Institutional Review Board of Sun Yat-sen University Cancer Center, and all methods were carried out in accordance with relevant guidelines and regulations or declaration of Helsinki. The need for informed consent was waived by the ethics committee of Sun Yat-sen University Cancer Center, because of the retrospective nature of the study.

Data availability statement

The data that support the findings of this study are available from the corresponding author ([email protected]) upon reasonable request.

Author contributions

B.D. contributed to conception and design; X.Z. contributed to development of methodology and contributed to acquisition of data; B.D., and X.Z., contributed to analysis and interpretation of data and writing of the manuscript; M.C. contributed to the study supervision and review of the manuscript. All authors approved the final version of the manuscript.

Siegel RL, Miller KD, Fuchs HE, Jemal A: Cancer Statistics, 2021. CA Cancer J Clin 2021, 71(1).
Brody H: Lung cancer. Nature 2020, 587(7834):S7.
Wu F, Wang L, Zhou C: Lung cancer in China: current and prospect. Curr Opin Oncol 2021, 33(1):40-46.
Miller M, Hanna N: Advances in systemic therapy for non-small cell lung cancer. BMJ 2021, 375:n2363.
Hanahan D, Weinberg RA: Hallmarks of cancer: the next generation. Cell 2011, 144(5):646-674.
Vander Heiden MG, Cantley LC, Thompson CB: Understanding the Warburg effect: the metabolic requirements of cell proliferation. Science 2009, 324(5930):1029-1033.
Hirschhaeuser F, Sattler UGA, Mueller-Klieser W: Lactate: a metabolic key player in cancer. Cancer Res 2011, 71(22):6921-6925.
Chang C-H, Qiu J, O'Sullivan D, Buck MD, Noguchi T, Curtis JD, Chen Q, Gindin M, Gubin MM, van der Windt GJW et al: Metabolic Competition in the Tumor Microenvironment Is a Driver of Cancer Progression. Cell 2015, 162(6):1229-1241.
Sun S, Li H, Chen J, Qian Q: Lactic Acid: No Longer an Inert and End-Product of Glycolysis. Physiology (Bethesda) 2017, 32(6):453-463.
Kumagai S, Koyama S, Itahashi K, Tanegashima T, Lin Y-T, Togashi Y, Kamada T, Irie T, Okumura G, Kono H et al: Lactic acid promotes PD-1 expression in regulatory T cells in highly glycolytic tumor microenvironments. Cancer Cell 2022, 40(2).
Wang Z-H, Peng W-B, Zhang P, Yang X-P, Zhou Q: Lactate in the tumour microenvironment: From immune modulation to therapy. EBioMedicine 2021, 73:103627.
Yoshida GJ: Metabolic reprogramming: the emerging concept and associated therapeutic strategies. J Exp Clin Cancer Res 2015, 34:111.
Gatenby RA, Gillies RJ: Why do cancers have high aerobic glycolysis? Nat Rev Cancer 2004, 4(11):891-899.
Tomida S, Takeuchi T, Shimada Y, Arima C, Matsuo K, Mitsudomi T, Yatabe Y, Takahashi T: Relapse-related molecular signature in lung adenocarcinomas identifies patients with dismal prognosis. J Clin Oncol 2009, 27(17):2793-2799.
Shedden K, Taylor JMG, Enkemann SA, Tsao M-S, Yeatman TJ, Gerald WL, Eschrich S, Jurisica I, Giordano TJ, Misek DE et al: Gene expression-based survival prediction in lung adenocarcinoma: a multi-site, blinded validation study. Nat Med 2008, 14(8):822-827.
Okayama H, Kohno T, Ishii Y, Shimada Y, Shiraishi K, Iwakawa R, Furuta K, Tsuta K, Shibata T, Yamamoto S et al: Identification of genes upregulated in ALK-positive and EGFR/KRAS/ALK-negative lung adenocarcinomas. Cancer Res 2012, 72(1):100-111.
Schabath MB, Welsh EA, Fulp WJ, Chen L, Teer JK, Thompson ZJ, Engel BE, Xie M, Berglund AE, Creelan BC et al: Differential association of STK11 and TP53 with KRAS mutation-associated gene expression, proliferation and immune surveillance in lung adenocarcinoma. Oncogene 2016, 35(24):3209-3216.
Pardoll DM: The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer 2012, 12(4):252-264.
Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, Li Z, Traugh N, Bu X, Li B et al: Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018, 24(10):1550-1558.
Barbieri I, Kouzarides T: Role of RNA modifications in cancer. Nat Rev Cancer 2020, 20(6):303-322.
Maeser D, Gruener RF, Huang RS: oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021, 22(6).
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, Li B, Liu XS: TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res 2020, 48(W1):W509-W514.
Ye F, Chen Y, Xia L, Lian J, Yang S: Aldolase A overexpression is associated with poor prognosis and promotes tumor progression by the epithelial-mesenchymal transition in colon cancer. Biochem Biophys Res Commun 2018, 497(2):639-645.
Fu H, Gao H, Qi X, Zhao L, Wu D, Bai Y, Li H, Liu X, Hu J, Shao S: Aldolase A promotes proliferation and G/S transition via the EGFR/MAPK pathway in non-small cell lung cancer. Cancer Commun (Lond) 2018, 38(1):18.
Buzzai M, Bauer DE, Jones RG, Deberardinis RJ, Hatzivassiliou G, Elstrom RL, Thompson CB: The glucose dependence of Akt-transformed cells can be reversed by pharmacologic activation of fatty acid beta-oxidation. Oncogene 2005, 24(26):4165-4173.
Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, Huang C, Li J, Dong X, Zhou Y et al: Integrated Proteogenomic Characterization of HBV-Related Hepatocellular Carcinoma. Cell 2019, 179(2).
Sarode P, Zheng X, Giotopoulou GA, Weigert A, Kuenne C, Günther S, Friedrich A, Gattenlöhner S, Stiewe T, Brüne B et al: Reprogramming of tumor-associated macrophages by targeting β-catenin/FOSL2/ARID5A signaling: A potential treatment of lung cancer. Sci Adv 2020, 6(23):eaaz6105.
Wan X, Guan S, Hou Y, Qin Y, Zeng H, Yang L, Qiao Y, Liu S, Li Q, Jin T et al: FOSL2 promotes VEGF-independent angiogenesis by transcriptionnally activating Wnt5a in breast cancer-associated fibroblasts. Theranostics 2021, 11(10):4975-4991.
Murarka S, Martín-Gago P, Schultz-Fademrecht C, Al Saabi A, Baumann M, Fansa EK, Ismail S, Nussbaumer P, Wittinghofer A, Waldmann H: Development of Pyridazinone Chemotypes Targeting the PDEδ Prenyl Binding Site. Chemistry 2017, 23(25):6083-6093.
Dharmaiah S, Bindu L, Tran TH, Gillette WK, Frank PH, Ghirlando R, Nissley DV, Esposito D, McCormick F, Stephen AG et al: Structural basis of recognition of farnesylated and methylated KRAS4b by PDEδ. Proc Natl Acad Sci U S A 2016, 113(44):E6766-E6775.
Zimmermann G, Papke B, Ismail S, Vartak N, Chandra A, Hoffmann M, Hahn SA, Triola G, Wittinghofer A, Bastiaens PIH et al: Small molecule inhibition of the KRAS-PDEδ interaction impairs oncogenic KRAS signalling. Nature 2013, 497(7451):638-642.
Papke B, Murarka S, Vogel HA, Martín-Gago P, Kovacevic M, Truxius DC, Fansa EK, Ismail S, Zimmermann G, Heinelt K et al: Identification of pyrazolopyridazinones as PDEδ inhibitors. Nat Commun 2016, 7:11360.
Canovas Nunes S, De Vita S, Anighoro A, Autelitano F, Beaumont E, Klingbeil P, McGuinness M, Duvert B, Harris C, Yang L et al: Validation of a small molecule inhibitor of PDE6D-RAS interaction with favorable anti-leukemic effects. Blood Cancer J 2022, 12(4):64.
Chang Y-T, Yeh Y-S, Ma C-J, Huang C-W, Tsai H-L, Huang M-Y, Cheng T-L, Wang J-Y: Optimization of a multigene biochip for detection of relapsed and early relapsed colorectal cancer. J Surg Res 2017, 220:427-437.
Michalik L, Desvergne B, Wahli W: Peroxisome-proliferator-activated receptors and cancers: complex stories. Nat Rev Cancer 2004, 4(1):61-70.
Yuan H, Lu J, Xiao J, Upadhyay G, Umans R, Kallakury B, Yin Y, Fant ME, Kopelovich L, Glazer RI: PPARδ induces estrogen receptor-positive mammary neoplasia through an inflammatory and metabolic phenotype linked to mTOR activation. Cancer Res 2013, 73(14):4349-4361.
Zuo X, Deguchi Y, Xu W, Liu Y, Li HS, Wei D, Tian R, Chen W, Xu M, Yang Y et al: PPARD and Interferon Gamma Promote Transformation of Gastric Progenitor Cells and Tumorigenesis in Mice. Gastroenterology 2019, 157(1):163-178.
Ali SR, Jordan M, Nagarajan P, Amit M: Nerve Density and Neuronal Biomarkers in Cancer. Cancers (Basel) 2022, 14(19).
Zuo X, Xu W, Xu M, Tian R, Moussalli MJ, Mao F, Zheng X, Wang J, Morris JS, Gagea M et al: Metastasis regulation by PPARD expression in cancer cells. JCI Insight 2017, 2(1):e91419.
Wang X, Qian C, Yang Y, Liu M-Y, Ke Y, Qian Z-M: Phosphorylated Rasal2 facilitates breast cancer progression. EBioMedicine 2019, 50:144-155.
Koh S-B, Ross K, Isakoff SJ, Melkonjan N, He L, Matissek KJ, Schultz A, Mayer EL, Traina TA, Carey LA et al: RASAL2 Confers Collateral MEK/EGFR Dependency in Chemoresistant Triple-Negative Breast Cancer. Clin Cancer Res 2021, 27(17):4883-4897.
Pan Y, Tong JHM, Lung RWM, Kang W, Kwan JSH, Chak WP, Tin KY, Chung LY, Wu F, Ng SSM et al: RASAL2 promotes tumor progression through LATS2/YAP1 axis of hippo signaling pathway in colorectal cancer. Mol Cancer 2018, 17(1):102.
Hui K, Yue Y, Wu S, Gu Y, Guan B, Wang X, Hsieh J-T, Chang LS, He D, Wu K: The expression and function of RASAL2 in renal cell carcinoma angiogenesis. Cell Death Dis 2018, 9(9):881.
Hui K, Gao Y, Huang J, Xu S, Wang B, Zeng J, Fan J, Wang X, Yue Y, Wu S et al: RASAL2, a RAS GTPase-activating protein, inhibits stemness and epithelial-mesenchymal transition via MAPK/SOX2 pathway in bladder cancer. Cell Death Dis 2017, 8(2):e2600.
Keenan TE, Tolaney SM: Role of Immunotherapy in Triple-Negative Breast Cancer. J Natl Compr Canc Netw 2020, 18(4):479-489.
Emens LA: Breast Cancer Immunotherapy: Facts and Hopes. Clin Cancer Res 2018, 24(3):511-520.
Mittendorf EA, Zhang H, Barrios CH, Saji S, Jung KH, Hegg R, Koehler A, Sohn J, Iwata H, Telli ML et al: Neoadjuvant atezolizumab in combination with sequential nab-paclitaxel and anthracycline-based chemotherapy versus placebo and chemotherapy in patients with early-stage triple-negative breast cancer (IMpassion031): a randomised, double-blind, phase 3 trial. Lancet 2020, 396(10257):1090-1100.
Schmid P, Cortes J, Pusztai L, McArthur H, Kümmel S, Bergh J, Denkert C, Park YH, Hui R, Harbeck N et al: Pembrolizumab for Early Triple-Negative Breast Cancer. N Engl J Med 2020, 382(9):810-821.
Pauken KE, Dougan M, Rose NR, Lichtman AH, Sharpe AH: Adverse Events Following Cancer Immunotherapy: Obstacles and Opportunities. Trends Immunol 2019, 40(6):511-523.
Cascone T, McKenzie JA, Mbofung RM, Punt S, Wang Z, Xu C, Williams LJ, Wang Z, Bristow CA, Carugo A et al: Increased Tumor Glycolysis Characterizes Immune Resistance to Adoptive T Cell Therapy. Cell Metab 2018, 27(5).
Bretz AC, Parnitzke U, Kronthaler K, Dreker T, Bartz R, Hermann F, Ammendola A, Wulff T, Hamm S: Domatinostat favors the immunotherapy response by modulating the tumor immune microenvironment (TIME). J Immunother Cancer 2019, 7(1):294.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of a novel glycolysis-related prognostic signature for predicting prognosis and tumor microenvironment of lung adenocarcinoma

Status:

Version 1

Abstract

Figures

1. Introduction

2. Materials And Methods

2.1. Data acquisition and preprocessing

2.2. Identification of differential expression and genetic alterations in GRGs

2.3. Identification of GRGs associated with survival and hierarchical clustering

2.4. Construction and validation of a prognostic model

2.5. Evaluation of the risk model of GRGs as an independent prognostic factor in LUAD

2.6. Enrichment analysis

2.7. Assessing the response to immunotherapy and drug sensitivity analysis

2.8. Statistical analysis

3. Results

3.1. Genetic variation of glycolysis-related prognostic genes in LUAD

3.2. Evaluation of glycolysis-related patterns based on 18 regulators

3.3. TME cell infiltration characteristics of glycolysis-related patterns

3.4. Construction of the glycolysis-related clusters

3.5. Construction of the glycolysis-related prognostic signature

3.6. Correlation of the glycolysis-related prognostic signature with the TME and immunotherapy response

4. Discussion

5. Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1