Development and Verification of the Amino Metabolism-Related and Immune-Associated Prognosis Signature in Gliomas

Aberrant reprogramming of metabolism has been considered a hallmark in various malignant tumors. The metabolic changes of amino acid not only have dramatic effects in cancer cells but also influence their immune-microenvironment in gliomas. However, the features of the amino acid metabolism-related and immune-associated gene set have not been systematically described. The expression level of mRNA was obtained from The Cancer Genome Atlas database and the Chinese Glioma Genome Atlas database, which were used as training set and validation set, respectively. Different bioinformatics and statistical methods were combined to construct a robust amino metabolism-related and immune-associated risk signature for distinguishing prognosis and clinical pathology features. Constructing the nomogram enhanced risk stratification and quantified risk assessment based on our gene model. Besides this, the biological mechanism related to the risk score was investigated by gene set enrichment analysis. Hub genes of risk signature were identified by the protein–protein interaction network. The amino acid metabolism-related and immune-associated gene signature recognized high-risk patients, defined as an independent risk factor for overall survival. The nomogram exhibited a high accuracy in predicting the overall survival rate for glioma patients. Furthermore, the high risk score hinted an immunosuppressive microenvironment and a lower sensitivity of immune checkpoint blockade therapy and also identified PSMC5 and PSMD3 as novel biomarkers in glioma. In conclusion, a novel amino acid metabolism-related and immune-associated risk signature for predicting prognosis in glioma has been constructed and identified as two potential novel biomarkers.


INTRODUCTION
Metabolic reprogramming is critical for maintaining the survival of cancer cells and defined as a hallmark of cancer, which might be the consequence of oncogenic mutations (1). Amino acid metabolism also emerges as an important role in the metabolic reprogramming of cancer cells because of its function in redox balance, energy regulation, biosynthesis support, and so on (2). Amino acids and their derivatives can not only regulate cancer cells but also modulate the surrounding microenvironment, which enhances the malignancy and immunosuppression of tumors (3), for instance, arginine derivations could change the chromatin structure to regulate gene expression, which promotes the proliferation of cancer cells (4). Kynurenine, which is the catabolic product of tryptophan, induces the invasion of cancer cells and the immunosuppression of a tumor microenvironment (5) by binding to transcription factor aryl hydrocarbon receptor (AHR) (6). Moreover, the activation of AHR hampers the performance of dendritic cells and T cells, which play a role in anti-tumor (7). The metabolism of amino acid is varied in tumors and plays a significant role not only in the biological process of tumor cells but also in the tumor microenvironment, particularly the modulation of the immune. All these indicate a better understanding of the metabolism of amino acids, which will offer potentially effective targets for cancer therapy (8).
In our study, we focus on glioma which is a group of highly heterogeneous neurocutaneous tumors, accounting for about 26% of all intracranial tumors, and is the most deadly primary malignant type of brain tumor in adults (9). Although combined therapy has been developed, including precise surgical resection, adjuvant radiotherapy, and temozolomide chemotherapy, the overall survival remains poor and has not been significantly improved. Furthermore, long-term survival is rare (10). Although immunotherapy has made a breakthrough progress in the treatment of a variety of solid tumors, the specific effect of immunotherapy in glioma is still not clear (11). It has been found that the expression level of immunosuppressive factors such as PD-L1 and IDO/TDO were dramatically elevated in gliomas. As is well known, PD-L1 could limit the function of effect T cells, and the metabolite mediated by IDO/TDO promotes the development of an immunosuppressive microenvironment (12,13). In addition, the upregulation of Treg cells could exhausting cytotoxic T cells to reduce the damage of tumor cells and enhance the immunosuppressive effects in the glioma microenvironment (14). However, how amino acid metabolism influences prognosis and the immune process in glioma progression needs further systematic research.
In our study, the amino acid metabolism-related and immune-associated risk signature was defined as an independent risk factor for the prognosis of glioma patients. The decision tree strongly verified the risk-dependent subgroups, and the nomogram showed an extremely high accuracy. In addition, a high risk score hinted an immunosuppressive microenvironment and lower sensitivity of ICB therapy, and PSMC5 and PSMD3 were identified as novel biomarkers in glioma. In summary, we demonstrated a novel amino acid metabolism-related and immune-associated risk signature for predicting prognosis in patients with glioma and identified two potential novel biomarkers.

MATERIALS AND METHODS
The workflow of our analysis is shown in Figure 1, and specific details are explained in the following sub-sections.

Data Preparation and Collection
The expression of mRNA and the clinical information of patients were collected from 698 patients in the cancer genome atlas database. Consistently, 413 samples were collected from 693 samples in the Chinese Glioma Genome Atlas, part B, as a validation set. Moreover, the mRNA data of normal brain tissue was collected from the Genotype-Tissue Expression Project (GTEx).

Determination of the Immune Status Through Single-Sample Gene Set Enrichment Analysis
Using the single-sample Gene Sets Enrichment Analysis (ssGSEA) algorithm based on the transcriptome profiling data and corresponding immunity-related gene sets retrieved from MSigDB (15,16) and using ESTIMATE algorithm, we analyzed the estimation of stromal and immune cells in tumor tissues (17), which has been developed to measure stromal level (stromal score), cyto-infiltration degree (immune score), and tumor purity.

Construction of Amino Metabolism-Related and Immune-Associated Signatures
The R package "WGCNA" was used to construct a scale-free coexpression network to verify a gene module that is mostly related to amino metabolism and immune in glioma. To explore the most robust genes, the LASSO regression model was performed (18). Furthermore, the risk scores were calculated by multiplying gene expression by the regression coefficient acquired upon Lasso regression. Based on the median risk score, all cases were divided into high-or low-risk groups. risk score in the prediction of 5-year overall survival (OS) by analyzing the receiver operating characteristic (ROC) curve (19). Next, a nomogram according to related prognostic factors was constructed to quantitatively predict the 1-, 2-, and 3-year survival rate in glioma patients. The CIBERSORT package in R was used to evaluate differences in the frequencies of 22 immune cell types (including Tregs and CD4+ cells) in glioma. The CIBERSORT is widely used in evaluating the type of immune cell in the microenvironment through estimating relative subsets of known RNA transcript (CIBERSORT) software (https:// cibersort.stanford.edu/) (20). The computational method was used in the low-and high-risk groups to explore the correlation of the TICs in different groups.

Gene Set Enrichment Analysis
In the Molecular Signatures Database, Hallmark and C7 gene sets were downloaded, which were used as the target gene sets to investigate the gene sets associated with risk score in the whole transcriptome of all glioma samples in the TCGA performed by the software GSEA-3.0. (NOM p < 0.05 and FDR q <0.05 were considered significant).

Functional Annotation for Genes of Interest and Construction of PPI
To explore the gene ontology (GO) of selected genes, R package cluster Profiler package was used to explore the functions among genes of interest, with a cutoff criterion of adjusted p <0.05. The GO annotation that contains the three sub-ontologiesbiological process, cellular component (CC), and molecular function-can identify the biological properties of genes and gene sets for all organisms (21). The Online tool Search Tool for the Retrieval of Interacting Genes (STRING) was used to predict protein-protein interactions (PPI) and construct a PPI network of selected genes. Using the STRING database, genes with a score of 0.4 were chosen to build a network model visualized by Cytoscape (v3.7.2) (22). In a co-expression network, Maximal Clique Centrality (MCC) algorithm was reported to be the most effective method of finding hub nodes (19). The MCC of each node was calculated by CytoHubba, a plugin in Cytoscape (23). In this study, the genes with the top 10 MCC values were considered hub genes.

Verification of the Expression Patterns and the Prognostic Values of Hub Genes
To explore the potential reliability of the hub genes, the expression level of each hub gene between cancer and normal tissue was plotted as a box plot graph. Based on the TCGA database, Kaplan-Meier univariate survival analysis was performed by using the survival package in R software to explore the relationship between overall survival and diseasefree survival with hub genes in patients. In the study, all the patients selected for survival analysis should be with complete clinical information. Consequently, based on the median expression value of hub genes, these samples were divided into two subgroups. The survival-related hub genes with log-rank p <0.05 were regarded as statistically significant.

Human Tissue Samples
Normal brain tissues were collected from patients who suffered from a serious brain injury.

Construction of Weighted Gene Co-expression Modules
First, combining the mRNA data in TCGA glioma and GTEx database, 10,550 differentially expressed genes between glioma and normal brain tissue were detected. Among these genes, 260 amino acid metabolism-related genes were ensured (Supplementary Figure S1A). The KEGG and GO analyses confirmed that these nods were mainly related to the biological function and pathway of amino acid metabolism. The PPI network showed a strong co-expressed correlation among the genes (Supplementary Figures S1B-D). Then, the samples in the training set were hierarchically clustered in the immunityhigh (immunity-H) or immunity-low (immunity-L) group by ssGSEA (Supplementary Figures S2A, B). The box plot of the fraction of immune cells in glioma tissues was significantly different among immunity-H and immunity-L groups (Supplementary Figure S2C). Consistently, the stromal scores, immune scores, and ESTIMATE scores of glioma samples in the immunity-H group remarkably increased compared with those in the immunity-L group (Supplementary Figures S2D, F). Meanwhile, the tumor purity in the immunity-H group was significantly lower than in the immunity-L group (Supplementary Figure S2G). Besides this, patients in the immunity-H group had a significantly poorer prognosis than those in the other groups (Supplementary Figure S2H).
To find the correlation between amino acid metabolism-related genes and immune infiltration in TCGA glioma, gene coexpression networks were constructed from the TCGA glioma datasets with the WGCNA package. Two modules in the TCGA glioma were recognized, and a different color was assigned for each module (Figure 2A). Then, we created a heat map of moduleimmune relationships to evaluate the association between each module and different immune scores (high and low). The results of the module-immune relationships showed that the gray module had the highest association with the immune-high group (pink module: r = 0.27, p < 0.001) in TCGA glioma ( Figure 2B). The module membership and gene significance were highly correlated in the gray module ( Figure 2C).

Identification of a 12-Gene Risk Signature Associated With Amino Acid Metabolism and Immune in Glioma
To identify the amino acid metabolism-related and immunoassociated risk signature, the univariate Cox regression analysis was used to select 30 genes in the training set, which were related to the prognosis of patients ( Figure 3A). Thereafter, the most relevant biomarkers for prognosis were identified through the LASSO Cox regression model, and overfitting was counteracted by 10-fold cross-validation. As a result, the group of 12 genes (PSMC5, GLUD1, DHTKD1, OGDH, PSMF1, PSMD3, PSMB8, PSMB9, PSMD5, PSMD12, PSMC1, and PSMD6) was extracted according to LASSO coefficients ( Figures 3B, C). The median risk score was defined as the cutoff value to divide the training set into two subgroups, including high-and low-risk groups, and a significant difference was found in both the molecular and clinical characteristics between these subgroups ( Figure 3D).
At the same time, there were significant differences according to the risk signature values of age-stratified and WHO-gradestratified clinical samples in both the TCGA and CGGA cohorts ( Figures 4A, B, E, F). The molecular pathological diagnosis of glioma has been put forward in clinical practice. IDH wild type and 1p19q non-codeletion gliomas were all the poor prognostic factors and had an inadequate response to traditional radiotherapy or chemotherapy of glioma patients (24). Such being the case, the distribution of the 12-gene signature was explored based on IDH status-stratified clinical samples ( Figures 4C, G) and 1p/19q codeletion status ( Figures 4D, H). Overall, these results indicated that the risk score based on the gene signature was significantly associated with clinical features.

Development of the Risk Score Signature and Assessment of the Predicting Capacity
Based on the groups with high and low risk scores, Kaplan-Meier analysis was performed, and it showed that patients with high risk scores had dramatically reduced overall survival compared with patients with a low risk score in both TCGA and CGGA datasets ( Figures 5A, B). Besides this, as far as the 1-, 3-, and 5year overall survival is in question, the values of the area under the curve (AUC) of the ROC curve for the TCGA glioma cohort were 0.875, 0.933, and 0.854. Consistently, concerning 1-, 3-, and 5-year overall survival, the values of AUC for the CGGA cohort were 0.641, 0.678, and 0.687, respectively ( Figures 5C, D). Furthermore, the plots were listed to show the distribution of gene expression, risk score, and survival status basing on the amino acid metabolism-and immune-related signature in TCGA and CGGA ( Figures 5E, F). To further explore the significance of our model in evaluating prognosis independently, we performed a univariate analysis as well as a multivariate analysis, which showed that the value of the risk score might be defined as an independent factor to evaluate the prognosis of glioma patients in both TCGA and CGGA ( Table 1).

Combination of the Risk Signature and Clinicopathological Features Improves Risk Stratification and Survival Prediction
To better enhance the risk stratification of prognosis, we constructed a decision tree through patients with different grades of glioma from TCGA. As a result, the difference in overall survival was observed in subgroups with different risk scores ( Figure 6A). Developing individualized treatment for individual glioma patients is necessary. Consistently assessing the potential risk and prognosis for individual glioma patients is also important. Consequently, we built a nomogram with risk score as well as clinical pathology features, including IDH mutation and 1p19q ( Figure 6D). Besides this, the calibration analysis was performed to elevate the accuracy of our nomogram. The results showed that the prediction line of the nomogram was extremely close to the ideal performance (45°dotted line) ( Figures 6B, C).

The Differences in Immunocyte Infiltration Degree and Enrichment Plots of Immune-Related Gene Sets From Gene Set Enrichment Analysis Between Highand Low-Risk TCGA Cohorts
Next, to explore whether our risk score partly assessed the immune status of the tumor microenvironment, the relationship of amino acid metabolism-and immune-related gene signature with the immunocyte infiltration degree was explored in gliomas. Interestingly, our results indicated that M2 (Cor = 0.31; p = 8.8e−6) and Tregs (Cor = 0.169; p = 0.0093) were obviously positively related to risk score ( Figures 7A, B). Furthermore, NK cells (Cor = -0.39; p =1.9e−08) and CD4+ T cells (Cor = -0.24; p = 0.00058) ( Figures 7C, D) showed a negative correlation with the risk score. Immunotherapy is increasingly becoming an important part of tumor therapy and can significantly improve the prognosis of cancer patients in a variety of solid tumors (25). Hence, we detected the expression of immune checkpoints in subgroups with a high or low risk score. According to our gene model, the expression level of PD-L1, PD-1, and CTLA-4 was lower in  glioma patients with a high risk score (P < 0.05) (Figures 7E, F). This result showed that the high-risk-score group may be less sensitive to immunotherapy. Furthermore, glioma with a high risk score was obviously enriched in the downregulation of the effect immunity pathway. We found that high risk score had a negative relationship with T cell migration ( Figure 7G). It is well known that enhancing T cell infiltration in glioma may increase the response rates to immunotherapy and increase survival. Consistently, high risk score had a negative relationship with tumor necrosis factor function ( Figure 7H). These results indicated that risk score could predict an immunosuppressive micro-environment.

Identification of Hub Genes From Risk Signature as Biomarkers in Glioma
The PPI network among the overlapped genes was established by using the STRING database and performing GO and KEGG (Supplementary Figure S3). MCC algorithm of CytoHubba plugin was used to select hub genes of the PPI network, and the hub genes are listed in Supplementary Table S1. Basing on the MCC scores, we selected the top 10 highest-scored genes from hub genes, including ODC1, OAZ2, PSMD2, PSMD12, PSMC1, PSMC5, PSMD3, PSME3, PSMD10, and PSMD5. The expression levels of these genes were verified according to the TCGA database. Kaplan-Meier plotter and the expression level of the top 10 genes were performed as shown in Supplementary Figure  S4 and Supplementary Figure S5. Then, we performed a multivariate Cox analysis to evaluate the prognostic value of these genes in gliomas ( Figure 8A and Supplementary Table S2).
In addition, we performed RT-PCR in our clinical samples which include six normal brain tissues, 24 WHO grade II, and 55 GBM samples. Consistently, we found that the expression level of PSMD3 was positive with the grade of glioma. Inversely, the expression level of PSMC5 was negative with the grade of glioma ( Figures 8B, D). Moreover, the results of the Kaplan-Meier analysis indicated that PSMD3 was significantly associated with worse overall survival of the glioma patients (P < 0.05) ( Figure 8C). Conversely, PSMC5 was significantly associated with better overall survival of the glioma patients (P < 0.05) ( Figure 8E). All the results confirmed that PSMD3 and PSMC5 can be identified as potential biomarkers in glioma patients.

DISCUSSION
The reprogramming of amino acid metabolism in gliomas has been reported to contribute to the malignancy biological process of glioma, including proliferation, migration, and so on. A previous study has constructed an amino acid-related risk signature for gliomas, which could predict the survival and clinical features of patients (26). However, more and more studies revealed that amino acid metabolism is not only caused by oncogene alterations but also changed the surrounding tumor microenvironment (27). In our study, we focused on the amino acid metabolism and immune status which was explored by ssGSEA and ESTIMATE and confirmed the differential expression of genes in gliomas. Then, WGCNA was performed to identify amino acid metabolism and immuno-related gene modules based on the data from TCGA, and an amino acid metabolism and immune signature was constructed by LASSO Cox regression model. Subsequently, the prognostic value of the gene signature was validated in CGGA cohorts. Our risk score system could distinguish high-risk patients, indicating that it can act as a confidential risk factor in the complex subgroups of patients. In addition, a decision tree has been constructed to enhance risk stratification on the basis of WHO grade and risk score, which showed that risk score could act as a major determinant. Moreover, the generation of the nomogram was used to quantify risk assessment and survival probability, which could show higher accuracy and discrimination in survival prediction compared with the traditional characteristics.
Furthermore, our risk score could also provide valuable information on immune cell infiltration in a tumor microenvironment and reflect the effects of immunotherapy. In the high-risk-score subgroup, the infiltration of immunosuppressive cells like Tregs is higher than in the lowrisk-score subgroup. Conversely, effect immune cells are decreased in the high-risk-score subgroup. Consistently, a recent study has shown that amino acid metabolism can regulate immune cells in cancer (3,28,29). Our gene model might provide a clue of how the microenvironment was influenced by amino acid metabolism. In addition, immune checkpoint therapy has shown a great potential for treatment in diverse solid tumors (30). However, the therapeutic efficacy has not lived up to expectations in gliomas, and the specific mechanisms for the problem still need more research. Different genomic subtypes or molecular profiles are the main challenges in the response to PD-1/PD-L1 checkpoint blockades (31). In addition, the amino acid derivatives could promote the immunosuppressive microenvironment and even affect the expression of immune checkpoints in glioma (32,33). Interestingly, the expression of PD-L1 and CTLA-4 in the high-risk group was significantly lower than in the low-risk group in our study. These results might indicate better efficacy and greater sensitivity of anti-PD1 therapy in low-risk glioma patients.
In our study, we identified two biomarkers by estimating amino acid and immune status in gliomas on the basis of the expression of mRNA. PSMD3, also known as P58 or RPN3, is one of the members in the proteasome subunit S3 family, which acts as the non-ATPase subunit of the 19S regulator lid (34). PSMD3 is widely expressed in most tissues and defined as an oncogene in various cancers. WBC and neutrophil counts are related to the expression of PSMD3 (35). Additionally, PSMD3 is also related to the glucose-related features of carbohydrates and fatty acids from the diet (36,37). Besides this, the higher level of PSMD3 mRNA predicts a worse prognosis of acute myeloid leukemia patients, and PSMD3 promotes the progression of chronic myeloid leukemia by stabilizing NF-kB (38,39). Consistently, PSMD3 is upregulated in breast cancer compared with normal tissue, and patients with a higher expression level of PSMD3 are related with worse survival. PSMD3 is strongly associated with the expression of HER2, which can stabilize HER2 from degradation (40). PSMC5 is defined as a 19S regulatory component, and it could identify and transform ubiquitin-labeled proteins into the form of degradation which can be mediated by 20S complex (41). Interestingly, PSMC5 directly regulates transcription, for instance, it can influence the activity class II trans-activator to regulate the transcription of MHC class II (42). Besides this, it can also recruit p53 to the promoter of p21 to upregulate its expression, which can decrease the DNA damage mediated by ultraviolet (43), yet the specific functions of both biomarkers in gliomas remain unclear. More research are needed for further study. Finally, several limitations of our work should be mentioned, namely: (1) although two biomarkers have been identified, the potential function of these genes still remains unclear and should be explored in a future study, and (2) tumor heterogeneity is one of the most important features in gliomas (44), which also means different microenvironment features exist among diverse tumor sites, yet all of the data and information are collected from public databases, which makes it impossible to detect the immune status in the same or diverse tumor regions. As a result, this gene signature should be better validated in well-designed, multicenter, prospective studies.

CONCLUSIONS
In summary, the construction and validation of 12 amino acid metabolism-and immune-related genes have been defined as a prognostic signature. This prognostic signature can predict the prognosis of patients and help to select individualized therapeutic strategy in clinical practice, which provides a comprehensive perspective for clarifying the underlying mechanisms that determine the prognosis for glioma. In addition, our risk score model is associated with the immune status of glioma patients, which may imply the potential effect of immuno-therapy. Besides this, we also have identified PSMC5 and PSMD3 as new biomarkers in glioma.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Institutional Ethics Committee of the Faculty of Medicine at Renmin Hospital of Wuhan University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
QC and FY designed the research. YX, LY, and RG carried out the work. PH, QS, and ST analyzed the data and wrote the paper. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the National Natural Science Foundation of China (no. 82072764).

ACKNOWLEDGMENTS
We gratefully acknowledge The Cancer Genome Atlas pilot project and Chinese Glioma Genome Atlas, which made the genomic data and clinical data of glioma available.