Glycolysis Define Two Prognostic Subgroups of Lung Adenocarcinoma With Different Mutation Characteristics and Immune Infiltration Signatures

Increasing studies have proved that malignant tumors are associated with energy metabolism. This study was aimed to explore biological variables that impact the prognosis of patients in the glycolysis-related subgroups of lung adenocarcinoma (LUAD). The mRNA expression profiling and mutation data in large LUAD samples were collected from the Cancer Genome Atlas (TCGA) database. Then, we identified the expression level and prognostic value of glycolysis-related genes, as well as the fractions of 22 immune cells in the tumor microenvironment. The differences between glycolysis activity, mutation, and immune infiltrates were discussed in these groups, respectively. Two hundred fifty-five glycolysis-related genes were identified from gene set enrichment analysis (GSEA), of which 43 genes had prognostic values (p < 0.05). Next, we constructed a glycolysis-related competing endogenous RNA (ceRNA) network which related to the survival of LUAD. Then, two subgroups of LUAD (clusters 1 and 2) were identified by applying unsupervised consensus clustering to 43 glycolysis-related genes. The survival analysis showed that the cluster 1 patients had a worse prognosis (p < 0.001), and upregulated differentially expressed genes (DEGs) are interestingly enriched in malignancy-related biological processes. The differences between the two subgroups are SPTA1, KEAP1, USH2A, and KRAS among top 10 mutated signatures, which may be the underlying mechanism of grouping. Combined high tumor mutational burden (TMB) with tumor subgroups preferably predicts the prognosis of LUAD patients. The CIBERSORT algorithm results revealed that low TMB samples were concerned with increased infiltration level of memory resting CD4+ T cell (p = 0.03), resting mast cells (p = 0.044), and neutrophils (p = 0.002) in cluster 1 and high TMB samples were concerned with increased infiltration level of memory B cells, plasma cells, CD4 memory-activated T cells, macrophages M1, and activated mast cells in cluster 2, while reduced infiltration of monocytes, resting dendritic cells, and resting mast cells was captured in cluster 2. In conclusion, significant different gene expression characteristics were pooled according to the two subgroups of LUAD. The combination of subgroups, TMB and tumor-infiltrating immune cell signature, might be a novel prognostic biomarker in LUAD.

Introduction

Lung adenocarcinoma (LUAD) is a highly fatal cancer of the respiratory system worldwide, and the 5-year survival rate of LUAD is less than 17% (Hirsch et al., 2017). It is reported that application of biomarkers may provide effective prognostic values in LUAD. For example, CAV1 and DCN participate in regulating LUAD progression (Yan et al., 2019). Despite progress in molecular researches, more prognostic biomarkers need to be further explored in LUAD.

Glycolysis is one of the well-studied pathways of glucose metabolism. Normal mammalian cells are inhibited by glycolysis under aerobic conditions, while malignant tumor cells are active in glycolysis even under a sufficient-oxygen environment. This metabolic feature of aerobic glycolysis is called the Warburg effect, which is manifested by high glucose uptake rate, active glycolysis, and high lactic acid content of metabolites (Allard et al., 1994). Previous studies have confirmed that increased levels of glycolysis are involved in malignancy progression, such as proliferation, invasion, and migration (Ganapathy-Kanniappan and Geschwind, 2013; Jin et al., 2017). It is reported that downregulation of Barx2 promotes aerobic glycolysis and predicts a poor prognosis in non-small cell lung carcinoma (NSCLC) (Chen et al., 2018). Given that, understanding the mechanisms associated with glycolysis could be a major breakthrough for finding potential prognostic targets. In this study, we aimed to explore new biomarker strategies of glycolysis-related tumor subgroups, which was associated with tumor-infiltrating immune cells, a tumor mutational burden (TMB) in LUAD.

The immune system plays key roles in the surveillance and elimination of tumor cells. The malignant phenotype of tumors is determined not only by the internal activity of tumor cells but also by the tumor-infiltrating immune cells in the tumor microenvironment (Swann and Smyth, 2007), which can promote or suppress the development and growth of tumors (Shiao et al., 2011). Previous studies have shown that the infiltration by immune cells in tumor tissues may have a prognostic value (Fridman et al., 2012). Therefore, the signature of tumor-infiltrating immune cells may be regarded as a potential predictive prognostic biomarker.

Tumor cell mutations could change the function and expression of proteins, resulting in the appearance of tumor-specific neoantigens. T-cells then recognize these neoantigens, causing an antitumor response. Increasing the activation of immune cells through immunotherapy may remove immune-mediated tumor cells and ultimately improve the patients’ prognosis. Researchers have identified that high levels of TMB may be more responsive to immunotherapy (Rooney et al., 2015; Fancello et al., 2019), and KRAS mutation can be regarded as a predictive biomarker in lung cancer (Hainsworth and Anthony Greco, 2016). Thus, further studies of mutation signature should be explored to identify prognostic biomarkers for LUAD.

One single prediction method always makes the results inaccurate. Therefore, in the present study, the combination of TMB and tumor-infiltrating immune cells to assess prognosis in different tumor subgroups of LUAD could provide more accurate and effective biomarkers.

Materials and Methods

Data Collection, Processing, and Validation

We firstly gathered RNA-sequencing data and prognostic information of LUAD from the GDC tool¹. Then, we downloaded somatic mutation data of all patients processed by VarScan software, and the “maftools” package and “pheatmap” package in R were used to further analyze the mutation data. Subsequently, the “ConsensusClusterPlus” package was performed to construct consensus analysis and PCA was applied to verify the accuracy of the classification. Eventually, the “limma” package was performed to identify differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed genes (DEGs) with |logFC > 1| and p < 0.05. The independent cohort GSE72094 was downloaded from the Gene Expression Omnibus (GEO) database² for data validation. The dataset contains both expression levels of related genes and clinical information such as age, gender, survival status, and survival time.

Gene Set Enrichment Analysis (GSEA)

Gene set enrichment analysis³ (Subramanian et al., 2005) was used to explore whether the selected gene sets revealed significant differences between the tumors tissues and normal tissues. We obtained glycolysis-related data (KEGG GLYCOLYSIS GLUCONEOGENESIS, BIOCARTA GLYCOLYSIS PATHWAY, HALLMARK GLYCOLYSIS, and REACTOME GLYCOLYSIS) from the MsigDB database⁴. p < 0.05 was the standard to assess whether to further investigate.

Construction of a Glycolysis-Related Competing Endogenous RNA (ceRNA) Network

The miRTarBase⁵ (Chou et al., 2018) and miRDB⁶ (Chen and Wang, 2020) databases were used to obtain miRNA-targeted glycolysis-related genes. StarBase⁷ (Li et al., 2014) was performed to predict relationships between lncRNAs and miRNAs.

Functional Enrichment Analysis of DEGs and Protein–Protein Interaction (PPI) Network Construction Across Two Clusters

To figure out the functional enrichment of glycolysis-related DEGs, the “GOplot” package and “ClusterProfiler” package were used to further analyze by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and p < 0.05 was considered as statistically significant. Then, a protein–protein interaction (PPI) network was constructed using the STRING database (Szklarczyk et al., 2015) and the cytoHubba plugin of Cytoscape software (Shannon et al., 2003).

Evaluation of Tumor-Infiltrating Immune Cells

CIBERSORT is a new algorithm used for calculating fractions of 22 immune cells subsets via RNA-seq or microarray data. The immune cells included seven types of T cells and three types of B cells, NK cells and various myeloid cells. We selected the samples with p < 0.05 to elevate the accuracy.

Validation of 43 Significant Glycolysis-Related DEGs Using Quantitative Reverse Transcription-Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from 16HBE and H1299 cells using the TRIzol reagent (Invitrogen, Carlsbad, CA, United States). After the purity and concentration of the total RNA were determined, the total RNA was reverse transcribed into cDNA using the PrimeScript RT reagent kit (Accurate Biology). The qRT-PCR was performed using the SYBR Green Premix Ex Taq II (Accurate Biology). The PCR conditions were set as follows: 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 30 s for each specific primer. Finally, the relative mRNA expression levels of 43 genes were calculated using the 2^–ΔΔCT method. The primer sequences are listed in Supplementary Data Sheet 2.

Statistical Analysis

Statistical analysis was utilized by R software (version 4.0.3) and GraphPad Prism (version 7.0). The differences of expression levels between two groups were determined by Student’s t-test. We used the “survival” packages to process the survival analysis, and survival curves were visualized by the Kaplan–Meier plotter⁸, which was examined by the log-rank test. The difference of infiltrating immune cells between the high-TMB group and the low-TMB group was determined by unpaired t-test, as well as the expression level of the glycolysis-related genes between cluster 1 and cluster 2. Differences were considered significant at p < 0.05.

Results

Glycolysis Is Associated With the Tumorigenesis of LUAD

We firstly obtained the RNA sequencing data and clinical data of LUAD from TCGA. Then, gene set enrichment analysis (GSEA) was used to further analyze the relationship between glycolysis and LUAD. The GSEA results showed that BIOCARTA GLYCOLYSIS PATHWAY, HALLMARK GLYCOLYSIS, and REACTOME GLYCOLYSIS were significantly enriched in LUAD, which illustrated that glycolysis is involved in the tumorigenesis of LUAD (Figure 1A). According to the GSEA results, 255 glycolysis-related genes were extracted and visualized by a heat map (Figure 1B). The information of 255 glycolysis-related genes is summarized in Supplementary Data Sheet 2. In addition, the flow diagram of this study is illustrated in Figure 2.

FIGURE 1

Figure 1. Relationship between glycolysis and tumorigenesis in LUAD. (A) Enrichment plots of three gene sets which were differentiated between in LUAD and normal tissues via GSEA, p < 0.05 was considered to be statistically significant. (B) The expression level of 255 glycolysis-related genes in LUAD and normal tissues. Red is up-regulated and blue is down-regulated. Genes are in rows; samples are in columns. LUAD, lung adenocarcinoma; GSEA, gene set enrichment analysis.

FIGURE 2

Figure 2. Flow chart of this study. TCGA, the Cancer Genome Atlas; GSEA, Gene set enrichment analysis; GEO, Gene Expression Omnibus; ceRNA, competing endogenous RNA.

Screening the Prognostic Value of Glycolysis-Related Genes in LUAD

To elucidate the relationship between the glycolysis-related genes and the prognosis of LUAD, we assessed the prognostic values of these 255 glycolysis-related genes using the Kaplan–Meier plotter database and found that only 43 genes were significantly associated with overall survival (all p < 0.05) (Table 1) (43 representative figures are shown in Supplementary Figure 1).

TABLE 1

Table 1. Differentially expressed 43 glycolysis-related genes in lung adenocarcinoma.

Consensus Clustering of Glycolysis-Related Genes Identified Two Clusters of LUAD

Glycolysis plays crucial roles in the biological processes of LUAD, so we next visualized the expression levels and correlation of 43 prognostic glycolysis-related genes (Figure 3). Then, we constructed glycolysis subgroups of LUAD with different biological characters based on 43 glycolysis-related genes. Unsupervised consensus clustering analysis indicated that k = 2 was the optimal number of clusters (Figures 4A–C) and principal component analysis (PCA) was used to verify the classification by 43 glycolysis-related genes (Figure 4D). Finally, a total of 534 patients were divided into cluster 1 and cluster 2, and the results showed that cluster 1 was associated with worse overall survival in LUAD patients compared with cluster 2 (p < 0.001) (Figure 4E). The specific TCGA patient ID of each cluster is displayed in Supplementary Data Sheet 2.

FIGURE 3

Figure 3. The landscape of 43 glycolysis-related genes in LUAD. (A) The expression levels of 43 glycolysis-related genes in LUAD. The red color means tumor tissues and the blue color means the non-tumor tissues. Genes are in rows; samples are in columns. (B) Vioplots of 43 glycolysis-related genes in LUAD, p < 0.05 was considered to be statistically significant. Blue is normal tissues and red is tumor tissues. (C) Spearman correlation analysis of 43 glycolysis-related genes in LUAD. The red color means positive correlation and the blue color means negative correlation. LUAD, lung adenocarcinoma.

FIGURE 4

Figure 4. Identification of consensus clusters of glycolysis-related genes in LUAD from the TCGA database. (A) Consensus clustering matrix for k = 2. (B) Consensus clustering CDF for k = 2–9. (C) Relative change in the area under the CDF curve for k = 2–9. (D) PCA of the 43 glycolysis-related genes in LUAD. (E) Survival curves of each cluster in LUAD. LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; CDF, cumulative distribution function; PCA, principal component analysis.

Validation of Expression of 43 Glycolysis-Related Genes, Consensus Clustering, and Prognostic Value Using an Independent Cohort

In order to verify the feasibility of grouping, the gene expression profile of GSE72094 was used for further analyses. The heatmap plot was visualized to further exhibit the distribution of 43 prognostic glycolysis-related genes (Figure 5A). As shown in Supplementary Figure 2, the expression of 42 of 43 prognostic glycolysis-related genes was significantly elevated and one gene was downregulated in LUAD tissues compared with adjacent normal lung tissues in GSE72094 (both p < 0.05). Meanwhile, qRT-PCR was also used to verify the expression of these genes. As expected, the results showed that 39 of the 43 genes showed the consistency with the above results, while four genes showed that there is no significance between 16HBE and H1299 (Supplementary Figure 3), which could be caused by the difference between tissue and cell lines. Furthermore, unsupervised consensus clustering analysis and PCA results showed that these two clusters are meaningful and can be verified using GSE72094 (Figures 5B–E). The survival analysis revealed that the joint of glycolysis and gene expression of 43 genes had a significant correlation with the prognosis of LUAD patients (Figure 5F).

FIGURE 5

Figure 5. Validation of the gene expression and grouping in GSE72094. (A) Heat map of 43 glycolysis-related genes between LUAD and normal tissue. (B) Consensus clustering matrix for k = 2. (C) Consensus clustering CDF for k = 2–9. (D) Relative change in area under the CDF curve for k = 2–9. (E) PCA of the 43 glycolysis-related genes in LUAD. (F) Survival curves of each cluster in LUAD. LUAD, lung adenocarcinoma; TCGA, The Cancer Genome Atlas; CDF, cumulative distribution function; PCA, principal component analysis.

Construction of a 43-Glycolysis-Related Signature

After screening 43 genes, we also constructed a ceRNA network to explore potential upstream mechanisms which may account for the worse overall survival in cluster 1. Using the limma R package according to the standards (p < 0.05, | logFC| > 1), a total of 186 DElncRNAs (Supplementary Data Sheet 2) and 150 DEmiRNAs (Supplementary Data Sheet 2) were identified in LUAD. Firstly, the miRNAs targeted by 43 glycolysis-related genes were predicted via miRTarBase and miRDB. After overlapping the results, a total of 33 DEmiRNAs (Supplementary Data Sheet 2) were identified from TCGA; the Kaplan–Meier plotter database displayed that 16 (Supplementary Figure 4 and Table 2) of 33 DEmiRNAs were significantly associated with overall survival (all p < 0.05). Then, 61 DElncRNAs (Supplementary Data Sheet 2) targeted by 16 DEmiRNAs were predicted via StarBase and 13 (Supplementary Figure 5 and Table 3) of 61 DElncRNAs have prognostic values (all p < 0.05). Eventually, we constructed a 43-glycolysis-related ceRNA network in LUAD, including 13 DElncRNAs, 16 DEmiRNAs, and 43 glycolysis-related genes (Figure 6).

TABLE 2

Table 2. Differentially expressed miRNA targeted by 43 glycolysis-related genes with prognostic value in lung adenocarcinoma.

TABLE 3

Table 3. Differentially expressed lncRNA with prognostic value in lung adenocarcinoma.

FIGURE 6

Figure 6. Construction of a ceRNA network. A flow diagram of a 43-glycolysis-related-gene ceRNA network in LUAD. ceRNA, competing endogenous RNA; LUAD, lung adenocarcinoma; DEmiRNAs, differentially expressed miRNAs; TCGA, The Cancer Genome Atlas; DElncRNAs, differentially expressed lncRNAs.

Identification and Functional Enrichment Analysis of the DEGs in Each Cluster

To further comprehend the relationship between glycolysis subgroups of LUAD and different patient prognoses, we identified the DEGs in each cluster from TCGA. A total of 384 DEGs (Supplementary Data Sheet 2), including 94 upregulated genes in cluster 1 and 290 upregulated genes in cluster 2 were screened by the limma R package on the basis of p < 0.05 and | log2(FC)| > 1. Those upregulated DEGs in cluster 1 are downregulated in cluster 2, which means that the rest of the unmentioned downregulated DEGs in cluster 1 are also elevated in cluster 2. In fact, from different perspectives, the DEGs are also different. For this reason, we only displayed and discussed the upregulated genes in each cluster for their different prognosis. To better realize the function of clustering, we further performed functional enrichment analysis. Upregulated DEGs of cluster 1 were mainly enriched in immune-related biological processes, such as complement activation, classical pathway, humoral immune response mediated by circulating immunoglobulin, and immunoglobulin-mediated immune response (Figure 7A). More specifically, the immune-related functions in cluster 2 mainly include cell chemotaxis, myeloid leukocyte migration, leukocyte chemotaxis, and neutrophil chemotaxis (Figure 7C). The results of KEGG enrichment analysis showed that cluster 1 and cluster 2 were mainly enriched in metabolism-related signaling pathways (Figures 7B,D). Moreover, the functional annotation of DEGs displayed that 10 GO terms were statistically significant (p < 0.05) (Figures 7E–H and Tables 4, 5). To identify the biological modules of DEGs in cluster 1, the PPI network was generated by using STRING and visualized by Cytoscape software. The cytoHubba plugin will sort genes by core degree through a variety of algorithms, and genes with high degrees tend to be hub genes. Then, the top 10 hub nodes were identified by the cytoHubba plugin, including HIST1H3B, HIST1H1B, HIST1H1D, HIST1H1E, HIST1H4C, HIST1H2BL, HIST1H2AH, HIST1H2BO, HIST1H3E, and HIST1H2AM (Figure 7I).

FIGURE 7

Figure 7. Functional enrichment analysis of DEGs in LUAD. (A,C) GO terms of DEGs in cluster 1 and cluster 2. (B,D) KEGG pathways enriched for the DEGs in cluster 1 and cluster 2. (E,F) The correlation between statistically top 10 DEGs and their GO terms in cluster 1 and cluster 2. (G,H) Functional analysis of the genes with higher expression using GO terms in cluster 1 and cluster 2. The outer circle means the expression (logFC) of DEGs in each enriched GO terms: Blue dots indicated the upregulated DEGs. The inner circle means the significance of GO terms (log10-adjusted p values). All p < 0.05 was considered to be statistically significant. (I) Hub genes selected via cytoHubba. GO, Gene Ontology; DEGs, differentially expressed genes; LUAD, lung adenocarcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

TABLE 4

Table 4. The top 10 gene ontology (GO) terms of differentially expressed genes in cluster 1.

TABLE 5

Table 5. The top 10 gene ontology (GO) terms of differentially expressed genes in cluster 2.

Identification of the Mutation Profile Features From Each Cluster

The simple nucleotide variation data of 561 LUAD samples were downloaded from TCGA and processed by VarScan software. The landscape of mutation data was visualized using the “maftools” package. Mutation profile features indicated that missense mutation was the most common type in Variant Classification, single-nucleotide polymorphism formed the nucleus of variant type, and C>A accounted for more components than other single-nucleotide variants in cluster 1 and cluster 2. The top 10 mutated signatures were visualized by a horizontal histogram, and the results showed that the mutated SPTA1 (29%) and KEAP1 (29%) are unique in cluster 1 among the top 10 mutated genes, as well as USH2A (27%) and KRAS (24%) in cluster 2, which may be the underlying mechanism of grouping (Figures 8A,B). Besides, the variant allele frequency (VAF) of LUAD is visualized in Figures 8C,D and a lollipop plot displayed these unique mutation points on the protein structure (Figures 8E,F). A waterfall plot displayed the mutation information in each sample, and in cluster 1, seven genes were mutated by >30%: TP53 (58%), TTN (56%), RYR2 (40%), MUC16 (38%), ZFHX4 (38%), CSMD3 (37%), and LRP1B (32%). In cluster 2, five genes were mutated by >30%: TP53 (44%), MUC16 (39%), TTN (38%), CSMD3 (34%), and RYR2 (32%) (Figures 9A,B). Besides, mutated oncogenic pathways were also visualized in each cluster (Figures 9C,D).

FIGURE 8

Figure 8. Summary of mutation profiling of LUAD in the MAF file. (A,B) Statistical calculations of mutation types in cluster 1 and cluster 2. (C,D) VAF of LUAD in cluster 1 and cluster 2. (E,F) Lollipop plot of special mutated genes in cluster 1 and cluster 2. MAF, minor allele frequency; LUAD, lung adenocarcinoma; VAF, variant allele frequency.

FIGURE 9

Figure 9. Visualization of the mutation profile signature in each cluster. (A,B) Waterfall plot of mutation profiles in cluster 1 and cluster 2 in LUAD samples. The comments below mean the mutation types. The bar chart on the right shows the distribution of mutant types in the top 30 genes. (C,D) The oncogenic pathways of mutant genes in cluster 1 and cluster 2. TMB, tumor mutation burden; LUAD, lung adenocarcinoma.

Identification of Immune Cell Infiltration Signatures of Each Cluster

The functional enrichment analysis showed that DEGs were involved in cell chemotaxis and immune-related biological processes, which lead us further to analyze the correlation of TMB with immune signatures in each cluster. Based on the CIBERSORT algorithm, we revealed the proportions of different immune cells in specific patients in each cluster. The percentage of 22 types of tumor-infiltrating immune cell in each cluster was visualized (Figures 10A,B). Subsequently, we indicated the differential abundance of immune cells in the low-TMB and high-TMB groups via violin plots. Memory resting CD4+ T cells (p = 0.03), resting mast cells (p = 0.044), and neutrophils (p = 0.002) showed higher infiltrating levels in the low-TMB group in cluster 1 (Figure 10C). Moreover, the infiltration levels of memory B cells (p = 0.044), plasma cells (p = 0.048), activated memory CD4+ T cells (p < 0.001), resting NK cells (p = 0.008), M1 macrophages (p < 0.001), and activated mast cells (p = 0.019) were higher in the high-TMB group in cluster 2, However, monocytes (p = 0.021), resting dendritic cells (p < 0.001), and resting mast cells (p < 0.001) showed a higher infiltrating level in the low-TMB group (Figure 10D).

FIGURE 10

Figure 10. Illustration of immune cell infiltrates features between low-TMB and high-TMB groups in each cluster. (A,B) The percent of 22 types of tumor-infiltrating immune cell in cluster 1 and cluster 2. (C,D) The difference of 22 types of immune cells between high and low TMB groups in cluster 1 and cluster 2. The higher TMB and lower TMB were shown in red and green. TMB, tumor mutation burden.

Discussion

Lung adenocarcinoma is one of the prevalent malignant tumors with increasing incidence, mortality, and poor prognosis (Nakamura and Saji, 2014; Hirsch et al., 2017). Meanwhile, LUAD displayed unique metabolic characteristics with high glycolytic activity, which is one of the most crucial hallmarks of cancer (Hanahan and Weinberg, 2011). It is reported that PPP1R14B-AS1 overexpression was associated with poor prognosis in LUAD (Yang et al., 2020). Inhibition of glycolysis could regulate the cell survival of LUAD (Farah et al., 2012). However, the relationship between glycolysis-related gene signatures and LUAD prognosis is still unclear so far. In this study, 43 differentially expressed glycolysis-related genes were related to the overall survival for LUAD patients. An appropriate recognition of tumor subgroups is necessary as it may affect patients’ prognosis and consideration of molecular detection. On the basis of 43 prognostic genes, we classified two subgroups of LUAD (cluster 1 and 2) by applying unsupervised consensus clustering, which had different clinical prognosis values. Therefore, we analyzed the possible mechanisms of different prognosis from multiple perspectives.

A subgroup-specific molecular study transforms the biological characteristics of LUAD into clinical prognostic subgroups to find potential biomarkers. Patients with different subgroups of LUAD had significantly different survival prognoses. In the present study, survival analysis results showed that cluster 1 was associated with worse overall survival for LUAD patients (p < 0.001). To explore the underlying upstream molecular mechanism which may account for the worse overall survival in cluster 1, we constructed a gene signature related to glycolysis, including 13 DElncRNAs, 16 DEmiRNAs, and 43 glycolysis-related genes and then verified the prognostic value of this gene signature in LUAD (all p < 0.05). In addition, the difference of gene expression level, immune infiltration signatures, and gene mutation infiltration signatures, and gene mutation information in cluster 1 was compared with that in cluster 2. The transformation of metabolic pathways is usually regulated by specific gene expression, and genetic signatures composed of multiple genes have been developed to predict the prognostic value of clinical cancer patients. According to the different clusters, we screened that 94 DEGs were enriched in immune-related biological processes in cluster 1. The top 10 hub genes were selected to help identify subgroups, and we found that histone modifications, one of epigenetic modifications, participated in the occurrence of cancer, which also could be regarded as a prognostic biomarker (Dworkin et al., 2009). Among these top 10 hub genes, HIST1H3B/C-K27M-mutated tumors exhibit a pro-angiogenic/hypoxic signature in diffuse intrinsic pontine gliomas (Castel et al., 2015). Copy number variations of HIST1H1B were connected with cellular development and proliferation in melanoma (Fidalgo et al., 2016). The expression of HIST1H2AH was higher in esophageal squamous cell carcinoma tissues than in adjacent non-tumorous tissues (Wang et al., 2019). Ten histone-coding genes here might regulate tumorigenesis to affect prognosis in cluster 1 for LUAD patients.

Compared with a single prediction model, the effect of analysis results combined with multiple perspectives can provide a more accurate prediction. Mutation was regarded as a therapy target to improve the prognosis of multiply tumors. Previous researches have reported that targeting SGLT2 may intercept LUAD progression at early stages (Scafoglio et al., 2018). SETD2-mutated LUAD patients exhibited poor recurrence-free survival (Kadara et al., 2017). In this study, we found that the mutated SPTA1 (29%) and KEAP1 (29%) are unique in cluster 1 among the top 10 mutated genes, as well as USH2A (27%) and KRAS (24%) in cluster 2. Furthermore, dysregulation of signaling pathways can also change in cancer metabolism to affect patients’ prognosis (Huang et al., 2016). For example, inactivation of the Hippo pathway is connected with the occurrence of various tumors (Panciera et al., 2017). Notch can promote glycolytic metabolism in T cell acute lymphoblastic leukemia (Palomero et al., 2007). Overexpressed SKA3 correlates with poor prognosis through the EGFR-PI3K-Akt axis in LUAD (Hu et al., 2020). Overexpression of CHAP2 may inhibit LUAD cell proliferation and correlate with high survival rate via the WNT signal pathway (Shang et al., 2019). In this study, the enrichment pathways of mutant genes are mainly WNT, PI3K, NOTCH, and Hippo signaling pathways in each cluster, which demonstrated that mutated genes contribute to different prognostic effects of LUAD through these pathways. However, the detailed mechanisms need further investigation.

The highly acidic microenvironment produced by tumor glycolysis may affect the infiltration of immune cells to varying degrees, eventually leading to immune escape and tumor progression (Gill et al., 2016; Cascone et al., 2018). Meanwhile, immune cell infiltration has a double feature to affect tumor progression, not only inhibiting the occurrence of tumor but also playing a pro-tumor role (Terlizzi et al., 2014). It is reported that lymphocyte infiltration has been connected with improved survival in NSCLC (Zeng et al., 2016). Mutated genes may generate neoantigens that increase lymphocyte infiltration in the tumor microenvironment. In this study, we found that the high TMB group has higher fractions of memory B cells (p = 0.044), plasma cells (p = 0.048), activated memory CD4+ T cells (p < 0.001), resting NK cells (p = 0.008), M1 macrophages (p < 0.001), and activated mast cells (p = 0.019) in cluster 2, which indicated that glycolysis-related cluster 2 with high TMB contributed to an immunomodulatory tumor microenvironment. These results indicated that TMB could affect the immune cell infiltration signatures and high TMB tends to cause the chemotaxis of immune cells in LUAD.

Conclusion

In conclusion, this study is the first to report a 43-gene prognostic signature related to glycolysis in LUAD. Furthermore, we constructed a complicated signature, among which the combination analysis of TMB and tumor-infiltrating immune cells was used to assess prognosis in different tumor subgroups of LUAD. Our study may provide a novel sight to realize the mechanisms of glycolysis and identify original gene targets for LUAD patients in the future.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

CH designed the research, analyzed the data, and wrote the manuscript. M-YZ and RL analyzed the data. T-TL and J-PL performed the data analysis and interpretation of the data. Y-QQ assisted with the study design and revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by grants from the Major Scientific and Technological Innovation Project of Shandong Province (2018CXGC1212), the CSCO-Qilu Cancer Research Fund (Y-Q201802-014), the Medical and Health Technology Innovation Plan of Jinan City (201805002), and the Technology Development Project (6010120081).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.645482/full#supplementary-material

Footnotes

1. ^ https://portal.gdc.cancer.gov/
2. ^ http://www.ncbi.nlm.nih.gov/geo/
3. ^ http://www.broadinstitute.org/gsea/index.jsp
4. ^ http://software.broadinstitute.org/gsea/msigdb/
5. ^ http://mirtarbase.mbc.nctu.edu.tw/
6. ^ http://www.mirdb.org/
7. ^ http://starbase.sysu.edu.cn/
8. ^ https://kmplot.com/analysis/

References

Allard, M. F., Schonekess, B. O., Henning, S. L., English, D. R., and Lopaschuk, G. D. (1994). Contribution of Oxidative-Metabolism and Glycolysis to Atp Production in Hypertrophied Hearts. Am.J. Physiol. 267, H742–H750.

Google Scholar

Cascone, T., McKenzie, J. A., Mbofung, R. M., Punt, S., Wang, Z., Xu, C., et al. (2018). Increased Tumor Glycolysis Characterizes Immune Resistance to Adoptive T Cell Therapy. Cell Metab. 27, 977–987.e4. doi: 10.1016/j.cmet.2018.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Castel, D., Philippe, C., Calmon, R., Le Dret, L., Truffaux, N., Boddaert, N., et al. (2015). Histone H3F3A and HIST1H3B K27M mutations define two subgroups of diffuse intrinsic pontine gliomas with different prognosis and phenotypes. Acta Neuropathol. 130, 815–827. doi: 10.1007/s00401-015-1478-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, H., Zhang, M., Zhang, W., Li, Y., Zhu, J., Zhang, X., et al. (2018). Downregulation of BarH-like homeobox 2 promotes cell proliferation, migration and aerobic glycolysis through Wnt/beta-catenin signaling, and predicts a poor prognosis in non-small cell lung carcinoma. Thorac Cancer 9, 390–399. doi: 10.1111/1759-7714.12593

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., and Wang, X. (2020). miRDB: an online database for prediction of functional microRNA targets. Nucleic Acids Res. 48, D127–D131. doi: 10.1093/nar/gkz757

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, C. H., Shrestha, S., Yang, C. D., Chang, N. W., Lin, Y. L., Liao, K. W., et al. (2018). miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 46, D296–D302. doi: 10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Dworkin, A. M., Huang, T. H., and Toland, A. E. (2009). Epigenetic alterations in the breast: implications for breast cancer detection, prognosis and treatment. Semin. Cancer Biol. 19, 165–171. doi: 10.1016/j.semcancer.2009.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Fancello, L., Gandini, S., Pelicci, P. G., and Mazzarella, L. (2019). Tumor mutational burden quantification from targeted gene panels: major advancements and challenges. J. Immunother. Cancer 7:183. doi: 10.1186/s40425-019-0647-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Farah, I. O., Lewis, V. L., Ayensu, W. K., and Cameron, J. A. (2012). Therapeutic implications of the warburg effect: role of oxalates and acetates on the differential survival of mrc-5 and a549 cell lines. Biomed. Sci. Instrum. 48, 119–125.

Google Scholar

Fidalgo, F., Rodrigues, T. C., Silva, A. G., Facure, L., de Sa, B. C., Duprat, J. P., et al. (2016). Role of rare germline copy number variation in melanoma-prone patients. Future Oncol. 12, 1345–1357. doi: 10.2217/fon.16.22

PubMed Abstract | CrossRef Full Text | Google Scholar

Fridman, W. H., Pages, F., Sautes-Fridman, C., and Galon, J. (2012). The immune contexture in human tumours: impact on clinical outcome. Nat. Rev. Cancer 12, 298–306. doi: 10.1038/nrc3245

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganapathy-Kanniappan, S., and Geschwind, J. F. H. (2013). Tumor glycolysis as a target for cancer therapy: progress and prospects. Mol. Cancer 12:152. doi: 10.1186/1476-4598-12-152

PubMed Abstract | CrossRef Full Text | Google Scholar

Gill, K. S., Fernandes, P., O’Donovan, T. R., McKenna, S. L., Doddakula, K. K., Power, D. G., et al. (2016). Glycolysis inhibition as a cancer treatment and its role in an anti-tumour immune response. Biochim. Biophys. Acta 1866, 87–105. doi: 10.1016/j.bbcan.2016.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hainsworth, J. D., and Anthony Greco, F. (2016). Lung Adenocarcinoma with Anaplastic Lymphoma Kinase (ALK) Rearrangement Presenting as Carcinoma of Unknown Primary Site: recognition and Treatment Implications. Drugs Real World Outcomes 3, 115–120. doi: 10.1007/s40801-016-0064-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanahan, D., and Weinberg, R. A. (2011). Hallmarks of cancer: the next generation. Cell 144, 646–674. doi: 10.1016/j.cell.2011.02.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirsch, F. R., Scagliotti, G. V., Mulshine, J. L., Kwon, R., Curran, W. J. Jr., Wu, Y. L., et al. (2017). Lung cancer: current therapies and new targeted treatments. Lancet 389, 299–311. doi: 10.1016/S0140-6736(16)30958-8

CrossRef Full Text | Google Scholar

Hu, D. D., Chen, H. L., Lou, L. M., Zhang, H., and Yang, G. L. (2020). SKA3 promotes lung adenocarcinoma metastasis through the EGFR-PI3K-Akt axis. Biosci. Rep. 40:BSR20194335. doi: 10.1042/BSR20194335

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, S., Chong, N., Lewis, N. E., Jia, W., Xie, G., and Garmire, L. X. (2016). Novel personalized pathway-based metabolomics models reveal key metabolic pathways for breast cancer diagnosis. Genome Med. 8:34. doi: 10.1186/s13073-016-0289-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, L., Chun, J., Pan, C., Alesi, G. N., Li, D., Magliocca, K. R., et al. (2017). Phosphorylation-mediated activation of LDHA promotes cancer cell invasion and tumour metastasis. Oncogene 36, 3797–3806. doi: 10.1038/onc.2017.6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadara, H., Choi, M., Zhang, J., Parra, E. R., Rodriguez-Canales, J., Gaffney, S. G., et al. (2017). Whole-exome sequencing and immune profiling of early-stage lung adenocarcinoma with fully annotated clinical follow-up. Ann. Oncol. 28, 75–82. doi: 10.1093/annonc/mdw436

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. doi: 10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, H., and Saji, H. (2014). Worldwide trend of increasing primary adenocarcinoma of the lung. Surg. Today 44, 1004–1012. doi: 10.1007/s00595-013-0636-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Palomero, T., Sulis, M. L., Cortina, M., Real, P. J., Barnes, K., Ciofani, M., et al. (2007). Mutational loss of PTEN induces resistance to NOTCH1 inhibition in T-cell leukemia. Nat. Med. 13, 1203–1210. doi: 10.1038/nm1636

PubMed Abstract | CrossRef Full Text | Google Scholar

Panciera, T., Azzolin, L., Cordenonsi, M., and Piccolo, S. (2017). Mechanobiology of YAP and TAZ in physiology and disease. Nat. Rev. Mol. Cell Biol. 18, 758–770. doi: 10.1038/nrm.2017.87

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell 160, 48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Scafoglio, C. R., Villegas, B., Abdelhady, G., Bailey, S. T., Liu, J., Shirali, A. S., et al. (2018). Sodium-glucose transporter 2 is a diagnostic and therapeutic target for early-stage lung adenocarcinoma. Sci. Transl. Med. 10:eaat5933. doi: 10.1126/scitranslmed.aat5933

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, J., Wang, Z., Chen, W., Yang, Z., Zheng, L., Wang, S., et al. (2019). Pseudogene CHIAP2 inhibits proliferation and invasion of lung adenocarcinoma cells by means of the WNT pathway. J. Cell. Physiol. 234, 13735–13746. doi: 10.1002/jcp.28053

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Shiao, S. L., Ganesan, A. P., Rugo, H. S., and Coussens, L. M. (2011). Immune microenvironments in solid tumors: new targets for therapy. Genes Dev. 25, 2559–2572. doi: 10.1101/gad.169029.111

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Swann, J. B., and Smyth, M. J. (2007). Immune surveillance of tumors. J. Clin. Invest. 117, 1137–1146. doi: 10.1172/JCI31405

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Terlizzi, M., Casolaro, V., Pinto, A., and Sorrentino, R. (2014). Inflammasome: cancer’s friend or foe? Pharmacol. Ther. 143, 24–33. doi: 10.1016/j.pharmthera.2014.02.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., An, S., Wang, D., Ji, H., Geng, M., Guo, X., et al. (2019). Quantitative Proteomics Identify the Possible Tumor Suppressive Role of Protease-Activated Receptor-4 in Esophageal Squamous Cell Carcinoma Cells. Pathol. Oncol. Res. 25, 937–943. doi: 10.1007/s12253-018-0395-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, Y. L., Xu, Z. J., Qian, L., Zeng, S. S., Zhou, Y. Y., Chen, X., et al. (2019). Identification of CAV1 and DCN as potential predictive biomarkers for lung adenocarcinoma. Am. J. Physiol. Lung Cell. Mol. Physiol. 316, L630–L643. doi: 10.1152/ajplung.00364.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Zhang, Y., Miao, L., Liao, W., and Liao, W. (2020). LncRNA PPP1R14B-AS1 Promotes Tumor Cell Proliferation and Migration via the Enhancement of Mitochondrial Respiration. Front. Genet. 11:557614. doi: 10.3389/fgene.2020.557614

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D.-Q., Yu, Y.-F., OuQ, -Y., Li, X.-Y., Zhong, R.-Z., Xie, C.-M., et al. (2016). Prognostic and predictive value of tumor-infiltrating lymphocytes for clinical therapeutic research in patients with non-small cell lung cancer. Oncotarget 7, 13765–13781. doi: 10.18632/oncotarget.7282

PubMed Abstract | CrossRef Full Text | Google Scholar