Comprehensive Characterization of a Novel E3-Related Gene Signature With Implications in Prognosis and Immunotherapy of Low-Grade Gliomas

Gliomas, a type of primary brain tumor, have emerged as a threat to global mortality due to their high heterogeneity and mortality. A low-grade glioma (LGG), although less aggressive compared with glioblastoma, still exhibits high recurrence and malignant progression. Ubiquitination is one of the most important posttranslational modifications that contribute to carcinogenesis and cancer recurrence. E3-related genes (E3RGs) play essential roles in the process of ubiquitination. Yet, the biological function and clinical significance of E3RGs in LGGs need further exploration. In this study, differentially expressed genes (DEGs) were screened by three differential expression analyses of LGG samples from The Cancer Genome Atlas (TCGA) database. DEGs with prognostic significance were selected by the univariate Cox regression analysis and log-rank statistical test. The LASSO-COX method was performed to identify an E3-related prognostic signature consisting of seven genes AURKA, PCGF2, MAP3K1, TRIM34, PRKN, TLE3, and TRIM17. The Chinese Glioma Genome Atlas (CGGA) dataset was used as the validation cohort. Kaplan–Meier survival analysis showed that LGG patients in the low-risk group had significantly higher overall survival time than those in the high-risk group in both TCGA and CGGA cohorts. Furthermore, multivariate Cox regression analysis revealed that the E3RG signature could be used as an independent prognostic factor. A nomogram based on the E3RG signature was then established and provided the prediction of the 1-, 3-, and 5-year survival probability of patients with LGGs. Moreover, DEGs were analyzed based on the risk signature, on which function analyses were performed. GO and KEGG analyses uncovered gene enrichment in extracellular matrix–related functions and immune-related biological processes in the high-risk group. GSEA revealed high enrichment in pathways that promote tumorigenesis and progression in the high-risk group. Furthermore, ESTIMATE algorithm analysis showed a significant difference in immune and stroma activity between high- and low-risk groups. Positive correlations between the risk signature and the tumor microenvironment immune cell infiltration and immune checkpoint molecules were also observed, implying that patients with the high-risk score may have better responses to immunotherapy. Overall, our findings might provide potential diagnostic and prognostic markers for LGG patients and offer meaningful insight for individualized treatment.


INTRODUCTION
Glioma is a primary type of tumor that occurs in the brain and spinal cord. The World Health Organization (WHO) classification system categorizes gliomas from grade I (lowest grade) through grade IV (highest grade) according to the malignant histological features (Wesseling and Capper, 2018). Low-grade gliomas (LGGs), which are less aggressive than glioblastoma multiforme (GBM), mainly originate from astrocytes and oligodendrocytes. Patients with LGGs are categorized as WHO grade II-III gliomas. The standard management of patients with LGG primarily involves surgical resection followed by adjuvant radiotherapy and chemotherapy (Lombardi et al., 2020). Even after patients receive these standard clinical interventions, the highly invasive nature and heterogeneity of LGGs can still lead to high rates of tumor recurrence and noteworthy malignant progression (Brat et al., 2015;Xia et al., 2018;Gittleman et al., 2020). Furthermore, the prognosis of LGGs varies diversely due to tumor heterogeneity.
LGG patients (mean age 41 years) are proposed to have survival averaging approximately 7 years, which is a significant sign of poor prognosis (Claus et al., 2015). Therefore, it is imperative to gain a more comprehensive understanding of the pathogenesis of LGG, identify effective and reliable biomarkers that could predict clinical outcomes, and formulate optimum therapeutic strategies.
Posttranslational modifications (PTMs) refer to covalent processing events of proteins which occur after synthesis and are normally mediated by diverse enzymes. Ubiquitination is a crucial posttranslational modification of a protein. It is an ATPdependent reversible process mediated by the ubiquitinproteasome system (UPS), including E1 ubiquitin-activating enzymes, E2 ubiquitin-conjugating enzymes, E3 ubiquitinprotein ligases, and deubiquitinating enzymes (DUBs) (Reyes-Turcu et al., 2009;Schulman and Harper, 2009;Buetow and Huang, 2016;Stewart et al., 2016). The dysregulation of UPS is largely involved in numerous biological functions, including cell cycle progression, cell proliferation, apoptosis, gene transcription, metastasis, transcriptional regulation, signaling, and inflammation (Deng et al., 2020). Accordingly, abnormal ubiquitination may contribute to various human pathologies such as neurodegeneration (Stieren et al., 2011;Popovic et al., 2014), autoimmune responses (Zangiabadi and Abdul-Sater, 2022), and oncogenic processes (Rape, 2018). In the UPS, E3 ubiquitin ligase serves as the essential part of the ubiquitination process owing to its substrate specificity (Zheng and Shabek, 2017). UPS is stringently regulated by E3 ligases that convey the specificity of ubiquitination. In particular, ubiquitin molecules are transferred from ubiquitinconjugating enzymes to specific substrates by E3 ubiquitin-protein ligases. The misregulation of UPS led by mutations in E3-related genes (E3RGs) is highly correlated with poor prognosis of cancers (Seeler and Dejean, 2017). Accumulating studies have demonstrated the tremendous contribution of E3-related proteins in glioma pathogenesis. For instance, MYH9-mediated ubiquitination of GSK-3β promotes malignant progression and chemoresistance in glioma (Que et al., 2022). The degradation of TUSC2 induced by NEDD4 facilitates glioblastoma progression (Rimkus et al., 2022). PARK2, frequently mutated in glioma, acts as a tumor suppressor by boosting ubiquitination-dependent degradation of β-catenin, which results in attenuation of Wnt signaling (Veeriah et al., 2010;Lin et al., 2015). Tripartite motifcontaining protein 11 (TRIM11), overexpressed in glioma, promotes proliferation, invasion, migration, and glial tumor growth via the induction of EGFR (Di et al., 2013). In addition, many E3 ligases have been reported to play vital roles in glioma carcinogenesis via the regulating PI3K/Akt pathway, such as SCF β−TrCP (Warfel et al., 2011), TRAF6 (Feng et al., 2014), and TRIM21 (Lee et al., 2017). Based on the significant function of E3 ligases in cancer pathogenesis, therapeutics targeting UPS have shown promising effects in clinical trials against cancers, such as PROTAC (proteolysis-targeting chimeric) (Qi et al., 2021). Two PROTAC strategies targeting CDK4 and/or CDK6 have been tested in glioma cells and are expected in clinical trials soon (Zhao and Burgess, 2019). Given the crucial roles of E3-related genes in glioma, the mechanism underlying the relationship between E3-related genes and the prognosis of LGG needs to be further addressed.
In the present study, a comprehensive analysis of E3RGs in LGG was conducted. Transcriptome data and clinical data of LGG samples were downloaded from The Cancer Genome Atlas (TCGA) database. Differentially expressed genes (DEGs) with prognostic significance were screened and identified. Subsequently, an E3-related prognostic signature was constructed, and the nomogram based on the risk signature was established to predict the prognosis of LGG. Meanwhile, enrichment analyses of risk-related differentially expressed genes and substrates of the risk signature genes were performed to disclose the underlying mechanism of LGG progression. Finally, the correlation between the risk signature and tumor immunity was analyzed. Our work provided an effective clinical tool for LGG prognosis prediction and preliminarily explored the biological functions and immune processes involved in the signature and the relative regulatory networks.

Datasets and E3-Related Gene Acquisition
The transcriptome sequencing data and corresponding clinical data of primary LGG were obtained from the TCGA database (https:// portal.gdc.cancer.gov/) and the Chinese Glioma Genome Atlas (CGGA) database, respectively, (http://www.cgga.org.cn/). The TCGA LGG cohort was selected as the training set, which included 451 tumor samples and five normal brain samples. The CGGA cohort (DataSet ID: mRNAseq_693) was chosen as the validation set, containing 240 primary LGG patients. The samples with incomplete clinical information and overall survival < 30 days had been excluded. Count data from TCGA and FPKM data from CGGA were applied for further analysis. The GSE68848 and GSE4290 datasets procured from the Gene Expression Omnibus database (GEO; https://www.ncbi.nlm.nih.gov/geo/) were utilized to validate the expression of the signature genes. The 630 E3RGs utilized in this study were collected from the ESBL database (https://esbl.nhlbi.nih.gov/Databases/KSBP2/Targets/Lists/E3ligases/) and iUUCD2.0 database (http://iuucd.biocuckoo.org/) (Supplementary Table S1).

Identification of Differentially Expressed E3-Related Genes
E3-related DEGs between LGG tissues and normal brain tissues were analyzed using the "limma," "edgeR," and "DESeq2" R packages with the cut-off criteria of |log 2 FC| ≥ 1 and p < 0.05 (Robinson et al., 2010;Love et al., 2014;Ritchie et al., 2015). The raw count data of the TCGA LGG cohort were used as the input for limma, edgeR, and DESeq. Volcano plots were generated to display DEG distribution from three algorithms mentioned earlier with the "tinyarray" R package. Venn diagrams were intersected to obtain the overlapping enriched terms also with the "tinyarray" R package.

Identification of E3-Related Differentially Expressed Genes With the Prognostic Value
The CPM of genes and clinical information were used for the subsequent analyses. The univariate Cox regression analysis and log-rank statistical test with the cut-off criteria of p < 0.05 for E3related DEGs were applied to select the genes with prognostic significance, using "survival" and "survminer" packages in R software. Venn diagrams were used to present the intersection of the enriched genes from Cox regression analysis and log-rank analysis.

Construction and Validation of the Prognostic Signature
To minimize the overfitting high-dimensional prognostic E3-related DEGs, least absolute shrinkage and selection operator (LASSO) regression analysis was performed with the "glmnet" R package (Friedman et al., 2010). Multivariate Cox regression analysis was conducted to construct prognostic models with the R package "My.stepwise". Hazard ratios (HRs) and 95% confidence intervals (CIs) were reported where applicable. The median risk score was calculated to categorize the LGG patients into high-risk and low-risk groups, based on the following formula: Risk score n i 1 Coef i × Expr i , (where Coef i is the coefficient and Expr i is the expression level of each intersected gene).
A Kaplan-Meier survival curve was used to determine the differences in overall survival using the R package "survival". Time-dependent receiver operating characteristic (ROC) analysis was executed to evaluate the prognostic accuracy of the seven-E3RG risk signature using R package "timeROC" (Blanche et al., 2013). The survival analysis result was presented by the R package "tinyarray" and "patchwork."

Development and Evaluation of the Nomogram
To assess whether the seven-E3RG prognostic risk signature can be utilized as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed using R package "survival" and "My.stepwise." The nomogram was constructed with independent prognostic parameters using the "rms" R package, and the calibration curves were utilized to reflect the accuracy of the nomogram.

Risk-Related Differentially Expressed Gene Analysis
LGG patients were divided into high-risk and low-risk groups according to the median risk score. The raw count data of the TCGA LGG cohort were used as the input data for limma, edgeR, and DESeq. "limma," "edgeR," and "DESeq2" R packages with the cut-off criteria of |log 2 FC| ≥ 1 and p < 0.05 were applied to two groups for DEG screening. Volcano plots were generated to display DEG distribution from three algorithms. Venn diagrams were intersected to obtain the overlapping enriched terms.

Functional Enrichment and GSEA
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed utilizing the "clusterProfiler" package to predict the biological function and related pathways (Yu et al., 2012). The top five enriched terms in the biological process (BP), cellular component (CC), and molecular function (MF) were visualized using "ggplot2" and "enrichplot" R packages. KEGG pathways with P. adjust<0.01 were chosen for presentation. Gene set enrichment analysis (GSEA) was conducted using GSEA v4.2.3 software with hallmark gene sets.

Immune Microenvironment Analysis
The tumor immune microenvironment (TIME) cell infiltration characterization was evaluated by the single-sample gene set enrichment analysis (ssGSEA) with the "GSVA" R package (Hänzelmann et al., 2013). An ESTIMATE algorithm was used to evaluate the immune and stromal activity in the LGG microenvironment with the "estimate" R package (Yoshihara et al., 2013). The correlation analysis was performed based on Pearson's correlation coefficient and presented using R packages "corrplot" and "circlize" (Gu et al., 2014).

Protein-Protein Interaction Network Analysis
A PPI network was constructed based on the substrates of E3RG signature with required interaction score > 0.4 using the STRING database (https://cn.string-db.org/) and visualized using Cytoscape (version 3.9.1) (Shannon et al., 2003;Szklarczyk et al., 2019). The top 10 hub genes in the PPI network were identified using the MCC algorithm with the CytoHubba plugin in Cytoscape.

Identification of Differentially Expressed E3-Related Genes With the Prognostic Value in Low-Grade Gliomas
The overall study workflow is presented in Figure 1. The TCGA cohort was used as the training set. The transcriptome and clinical data from TCGA included 451 tumor samples and five normal brain samples. The 630 E3RGs were selected and applied in this study (Supplementary Table S1). Three differential expression analyses were performed. According to the | log 2 FC | > 1.0 and p < 0.05, DEGs were displayed in volcano plots ( Figure 2A). The overlapping of DEGs from three differential expression analyses indicated that 44 genes were upregulated and 47 genes were downregulated in tumor samples ( Figure 2B). Subsequently, DEGs were further analyzed by univariate Cox regression analysis and log-rank statistical test to evaluate the prognostic significance. In total, 38 E3RGs with a significant prognostic value were obtained by the intersection of results from both analyses ( Figure 2C).

Construction of E3-Related Gene Prognostic Risk Signature in Low-Grade Gliomas
To further explore the prognostic value of E3RGs in LGGs, 38 overall survival-associated E3RGs were incorporated into the LASSO regression (Liu et al., 2021;2022a;2022b;2022c), 20 of which were selected for further multivariate Cox regression analysis ( Figures 3A,B). Following this, a seven-E3RG prognostic signature significantly correlated with LGG prognosis was developed by performing multivariate Cox regression analysis and shown in the forest plot ( Figure 3C). The risk score was calculated for each LGG patient by the following formula: The expression of the signature genes was validated in both GSE4290 and GSE68848 datasets and proved consistent with that in the TCGA cohort ( Figures 3D,E). Univariate Cox regression analysis revealed that seven signature genes were all strongly correlated with the prognosis ( Figure 4A). As shown in the Kaplan-Meier curves, four genes of the E3RG signature were considered to have favorable prognostic effects, while three were considered to have poor prognostic effects ( Figure 4B).
Based on the seven-E3RG risk signature, patients were divided into the high-risk and low-risk subgroups according to the median risk score in TCGA cohorts. The resulting Kaplan-Meier curve displayed a significant difference in overall survival between the LGG patients in the high-risk and the low-risk group, suggesting the established signature effectively predicts survival ( Figure 5A; p < 0.0001). The overall survival of patients with low-risk scores was significantly higher than that of patients with high-risk scores. The genes referred to in the signature were remarkably differentially expressed in the high-risk group and the low-risk group. The

Validation of the Risk Signature in the Chinese Glioma Genome Atlas Cohort
To test whether the prognostic gene signature has similar predictive performance and accuracy in other LGG cohorts, the CGGA cohort was used as a validation set. In the CGGA cohort, the patients were divided into low-risk and high-risk groups by the median risk score with the same formula calculated from the TCGA cohort ( Figure 6B). Consistent with the results obtained from the training set, survival analysis using the Kaplan-Meier method exhibited a better prognosis for patients in the low-risk group ( Figure 6A, p < 0.0001). The distribution of risk score and survival status is shown in Figures 6B,C. A heatmap of gene expression in the CGGA cohort is presented in Figure 6D, based on the risk score. The predicted AUCs of 1 year, 3 year, and 5 year are 0.80, 0.79, and 0.71, respectively ( Figure 6E). Prognostic analyses showed similar results. These results demonstrated that the E3RG risk signature was positively correlated with LGG prognosis.

Independent Prognostic Value of the E3-Related Gene Risk Signature in the Cancer Genome Atlas and the Chinese Glioma Genome Atlas Low-Grade Glioma Cohorts and Construction of a Nomogram
We first evaluated the prognostic value of age, gender, and risk score in patients with LGG from the TCGA cohort through univariate Cox regression analysis. The result revealed that both age (HR = 1.07, 95% CI = 1.05-1.09, p < 0.001) and risk score (HR = 2.72, 95% CI = 2.26-3.27, p < 0.001) of LGG patients were significantly correlated with overall survival ( Figure 7A). Moreover, multivariate Cox regression analysis showed that age (HR = 1.05, 95% CI = 1.03-1.07, p < 0.001) and risk score (HR = 2.28, 95% CI = 1.87-2.76, p < 0.001) affected overall survival as independent prognostic factors ( Figure 7B). Similar results were obtained when univariate and multivariate Cox regression analyses were applied in LGG patients from the CGGA cohort. Of note, gender (HR = 1.56, 95% CI = 1-2.44, p = 0.049) and risk score (HR = 1.97, 95% CI = 1.62-2.38, p < 0.001) were correlated with overall survival, but only risk score (HR = 1.97, 95% CI = 1.62-2.39, p < 0.001) became an LGG samples and five normal brain samples identified using edgeR, limma, and DESeq2 algorithms, with the cut-off criterion p < 0.05 and |log 2 FC| ≥ 1. Blue dots: significantly downregulated genes; red dots: significantly upregulated genes. (B) Venn diagram of the overlapping E3-related DEGs screened by the three differential expression analyses. (C) Venn diagram of the intersection for overall survival-correlated E3-related DEGs identified using univariate Cox regression analysis and log-rank statistical test with the cut-off criteria of p < 0.05. Left, predicted favorable prognosis; right, predicted poor prognosis.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 905047 independent prognostic factor in the CGGA cohort ( Figures 7C,D). The subtle difference may be attributed to the ethnicity difference.
On the basis of the seven-E3RG risk signature, we established a nomogram that could predict the prognosis of LGG in the TCGA cohort ( Figure 7E). Briefly, the points of different variables were mapped to the corresponding lines, while the total points of the patients were calculated and normalized to a distribution of 0-100. By performing this, the 1-, 3-, and 5-year survival status for LGG patients could be approximately estimated based on the prognosis axis and total point axis. In addition, the calibration curves for the probability of 1-, 3-, and 5-year overall survival showed a strong consistency between the predicted value of the nomogram and the actual value in both the TCGA and CGGA cohorts ( Figures 7F,G). Thus, the nomogram could serve as a favorable reference for clinical decision-making.

Identification and Function Analysis of Risk-Related Differentially Expressed Genes
To further investigate the potential biological functions and pathways of the risk signature, we screened the DEGs by three Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 905047 differential expression analyses between the high-risk group and the low-risk group in the TCGA LGG cohort. These results were displayed in volcano plots ( Figure 8A). As shown by the Venn diagrams in Figure 8B, 528 upregulated genes and 134 downregulated genes in the high-risk group were identified and applied for GO and KEGG pathway analyses ( Figure 8B). From the GO analysis, the risk-related DEGs were enriched in extracellular matrix-related functions, which suggested stronger migration and invasion potentials of tumors for LGG patients in the high-risk group, such as collagen-containing extracellular matrix and extracellular matrix structural constituent. Meanwhile, the risk-related DEG high enrichment was observed in several immune-related biological processes, such as MHC class II protein complex, MHC protein complex, MHC class II receptor activity, and immune receptor activity ( Figure 8C). Moreover, the KEGG pathway enrichment analysis showed that risk-related DEGs were principally intensified in immune-related pathways and cell adhesion molecules, which reflected the malignant characteristics of LGG in the high-risk group ( Figure 8D). This was also consistent with the results of GO analysis.

GSEA Enrichment Analysis
To elucidate the potential functional differences between the high-risk and low-risk groups, GSEA was performed with the TCGA LGG cohort. GSEA revealed that pathways related to inflammatory response, such as complement, IL-2/STAT5 signaling, IL-6/JAK/STAT3 signaling, and inflammatory response, were enriched in the high-risk group. Pathways that promote tumorigenesis and progression were also enriched in the high-risk group, such as PI3K/AKT/MTOR signaling, mTORC1 signaling, epithelial-mesenchymal transition, glycolysis, and KRAS signaling up. The enrichment of genes related to E2F targets, G2M checkpoint, and mitotic spindle suggested the correlation of cell cycle dysregulation with the risk score ( Figure 9).

The Role of the E3-Related Gene Risk Signature in the Tumor Immune Microenvironment
In order to further investigate the roles of the risk signature in TIME cell infiltration, we evaluated the landscape of 28 TIMEinfiltrating cell types in the low-risk and high-risk groups by ssGSEA (Charoentong et al., 2017). In total, 25 TIME cell types presented significant differences in infiltration between low-risk and high-risk groups ( Figure 10A). Although eosinophils, monocytes, and CD56 dim natural killer cells were not highly enriched in the high-risk group, a mild increase was still noted in eosinophils and monocytes. The expression of immune cell markers was displayed in a heatmap ( Figure 10B). Next, the ESTIMATE algorithm was applied to evaluate the immune and stromal activity in the LGG tumor microenvironment. The results disclosed that the immune and stromal activities were significantly elevated in the high-risk group, which might provide evidence for the contribution of an inflammatory environment of LGG as well ( Figures 10C,D). Correlation analysis was performed, and potential relations between the risk score signature and each TIME cell type are shown in Figure 10E. A significant positive correlation between the risk score and TIME cell infiltration was observed, except for eosinophils, monocytes, and CD56 dim natural killer cells. The risk score was proven to be positively correlated with the expression of immune checkpoint molecules ( Figure 10F), implying the meaningful roles of risk score signature in predicting the possible response of LGG patients to clinical immunotherapy.

Function Analysis of Substrates of E3-Related Gene Risk Signature in Low-Grade Gliomas
To acquire a better understanding of the potential biological function of the risk signature, the potential substrates of E3related gene signature were searched on Ubibrowser 2.0 (http:// ubibrowser.bio-it.cn/ubibrowser_v3/). The known substrates were applied for protein-protein interaction network analysis using the STRING database (https://cn.string-db.org/). A PPI network of 76 substrates of the risk signature, including 69 nodes and 772 edges, was constructed using the STRING database ( Figure 11A). The top 10 hub genes with the highest linkage degrees were obtained using the MCC algorithm of the cytoHubba plugin in Cytoscape3.9.1. These genes included PARK2, VDAC1, DNM1L, MFN2, MFN1, PINK1, TOMM20, PARK7, BCL2L1, and SNCA ( Figure 11B). To disclose the potential biological functions of the substrates that were involved, GO and KEGG pathway analyses were performed. GO analysis presented high enrichment in neuron death (regulation of neuron death, neuron death, negative regulation of neuron death, apoptotic mitochondrial changes, and death domain binding) and ubiquitin-related functions (ubiquitin-like protein ligase binding and ubiquitin-protein ligase binding) ( Figure 11C). KEGG pathway analysis results demonstrated that the substrates were mostly correlated with neurodegeneration (pathways of neurodegeneration, Parkinson's disease, and amyotrophic lateral sclerosis), cell death (mitophagy, apoptosis, and autophagy), and immune response (NOD-like receptor signaling pathway, PD-L1 expression, and PD-1 checkpoint pathway in cancer) ( Figure 11D).

DISCUSSION
In the past few decades, overwhelming evidence indicates that E3 ubiquitin ligases play pivotal roles in tumorigenesis, cancer progression, and treatment responses. Since E3 ligases determine the targets of the UPS, they play an essential role in cellular functions. They take part in biological processes, including but not limited to apoptosis, cell growth, senescence, proliferation, immune system evasion, metabolism, DNA repair, inflammation, invasion, metastasis, and angiogenesis. In GBM, alterations in EGFR are commonly seen (Brennan et al., 2013). EGFR could be downregulated by PARK2 but increased by TRIM11; both are E3 ligases possessing opposite effects on EGFR (Di et al., 2013;Lin et al., 2015). The PI3K/Akt pathway is altered in approximately 90% of GBM cases (Brennan et al., 2013). The PI3K/Akt signaling could be modulated by the SCF β−TrCP complex and SCF Skp2 complex and regulate the proliferation of primary GBM cell lines, glioma stem cells (GSCs), and established GBM cell line models (Winston et al., 1999;Li et al., 2009;Yang et al., 2009;Chan et al., 2012;Feng et al., 2014). Hence, more and more E3 ubiquitin ligases emerge as potential targets of drug designs for cancer therapies as they own better specificity for the recognition of substrates.
In this study, we first identified 38 DEGs with survival significance, in which seven DEGs showed strong prognostic performance and constituted a risk signature. The signature consists of AURKA, MAP3K1, TRIM34, PCGF2, PRKN, TLE3, and TRIM17. Patients in the high-risk group were more likely to have a worse prognosis compared with the ones in the low-risk group. Among the seven genes of the risk signature, many have been reported in glioma pathogenesis. Aurora kinase A (AURKA) has emerged as a drug target for glioblastoma for being highly involved in cell proliferation, migration, and invasion (Wang et al., 2018;Nguyen et al., 2021). MAP3K1 might promote glioma stem cell progression and be positively associated with resistance to temozolomide (TMZ) and radiotherapy (Bi et al., 2020;Wang et al., 2020). PRKN, first found to be mutated in patients with early-onset Parkinson's disease, has also been confirmed to carry mutations and deletions in human malignancies including glioblastoma, colon cancer, and lung cancers (Veeriah et al., 2010). In addition, PRKN inhibits glioma cell growth in vitro and in vivo by downregulating the intracellular levels of β-catenin and EGFR, leading to decreased activation of both Wnt-and EGFstimulated pathways (Lin et al., 2015). Although having not been reported in glioma-related research, TRIM34, PCGF2, TLE3, and TRIM17 have been revealed to contribute to the carcinogenesis of other cancers. TRIM34 appears to attenuate colon inflammation and tumorigenesis by sustaining barrier integrity, highlighting its role in immune responses (Lian et al., 2021). PCGF2 serves as a tumor suppressor in breast cancer, gastric cancer, and colon cancer probably for the negative regulation of Akt activation Guo et al., 2010;Zhang et al., 2010). TLE3 expression is positively correlated with taxane sensitivity in patients with ovarian carcinoma but not breast cancer (Samimi et al., 2012;Bartlett et al., 2015). TRIM17 augments BRAF-targeted therapy sensitivity of melanoma cells by preventing BCL2A1 from being ubiquitinated and degraded by TRIM28 (Lionnard et al., 2019). The studies mentioned earlier once again We screened the DEGs between high-risk and low-risk groups and performed the GO and KEGG pathway analyses. GO analysis revealed that the DEGs were enriched in extracellular matrix-related (ECM-related) functions in the high-risk group. The glioma ECM has several unique characteristics that make it distinct from the ECM of normal brain tissue. Glioma cells express components such as tenascin-C, fibronectin, and thrombospondin, which support the adhesion and migration of glioma cells. Furthermore, signals driven by the ECM also help shape tumor phenotypes (Sood et al., 2019). Our result corroborated that glioma ECM is tightly related to tumor cell proliferation and differentiation and poor prognosis. In addition, DEG high enrichment was observed in several immune-related biological processes such as MHC-related complex and immune receptor activity. The glioma TIME has diverse cell types and immune cell infiltration which consequently create a field of dynamic cytokine and chemokine communication. MHCII is expressed in different types of gliomas and is associated with increased infiltration of FIGURE 7 | regarding overall survival in the TCGA cohort (A,B) and the CGGA cohort (C,D) (*p < 0.05; ***p < 0.001). (E) Construction of a nomogram based on the independent prognostic values in the TCGA cohort. The calibration curves between predicted and observed 1-year, 3-year, and 5-year outcomes of nomograms in the TCGA cohort (F) and the CGGA cohort (G). Gray diagonal line represented ideal prediction.
FIGURE 8 | Identification and functional enrichment analyses of risk-related DEGs. (A) Volcano plot of risk-related DEGs between the high-risk group and low-risk group identified using edgeR, limma, and DESeq2 algorithms, with the cut-off criterion p < 0.05 and |log 2 FC| ≥ 1. Blue dots: significantly downregulated genes; red dots: significantly upregulated genes. (B) Venn diagram of the overlapping risk-related DEGs screened by the three differential expression analyses. (C) GO analysis of riskrelated DEGs with three terms biological processes, cellular components, and molecular functions (P.adjust < 0.05). (D) KEGG pathways enriched in the upregulated and downregulated risk-related DEGs (P.adjust < 0.01).
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 905047 T cells. Both clinical and transcriptomic data have uncovered that tumoral MHCII is tightly correlated with poor prognosis and immune responses (Chih et al., 2021, 01). The KEGG pathway enrichment analysis showed that risk-related DEGs were primarily intensified in immune-related pathways and cell adhesion molecules. Accordingly, both pathways contribute to the invasion and metastasis of glioma. GSEA results showed that pathways related to inflammatory responses were enriched in the high-risk group. Previous research indicates that about 15-20% of cancer cases suffer infection, chronic inflammation, or autoimmunity in the same tissue before solid tumor formation (Grivennikov et al., 2010). Meanwhile, the inflammatory nature of the tumor microenvironment promotes the development and survival of tumors (Mantovani et al., 2008). IL-2/STAT5 signaling is crucial for the regulation of regulatory T cells and immune tolerance (Shi et al., 2018). Similarly, the IL6-JAK-STAT3 signaling pathway may promote tumor cell proliferation, invasion, and metastasis and suppress immune response (Johnson et al., 2018). In addition, inflammatory responses and complement-related pathways could affect the tumorigenesis, immune surveillance, and immunotherapy response (Greten and Grivennikov, 2019). Our study showed the enrichment of the immune-related pathways in the high-risk group, which emphasized their significant roles in glioma prognosis. Furthermore, pathways that facilitate tumorigenesis and progression were also enriched in the high-risk group. The PI3K/AKT/mTOR pathway is widely dysregulated almost in all human cancers and is pivotal to cancer cell proliferation, survival, and therapy resistance (Cirone, 2021;Pungsrinont et al., 2021). Mutations within genes of this signaling pathway are the most common events occurring in solid malignancies including glioma (Baghery Saghchy Khorasani et al., 2021). mTORC1, one of the mTOR forms, comprises mTOR, raptor, GβL, and deptor. In particular, mTORC1 signaling is mainly involved in cell growth and metabolism (Unni and Arteaga, 2019). Furthermore, AKT/PI3K signaling could indirectly activate mTORC1 by the phosphorylation of PRAS40, a known mTORC1 inhibitor (Sancak et al., 2007;Thedieck et al., 2007;Vander Haar et al., 2007;Wang et al., 2007). Activated mTORC1 signaling may trigger recurrent reprograming that helps escape from glycolytic addiction in cancer cell lines from various solid tumor types (Pusapati et al., Recently, treatments targeting PI3K/AKT/mTOR and mTORC1 pathways have emerged as promising strategies in cancer therapeutics (Yang et al., 2019;Peng et al., 2022). Epithelial-mesenchymal transition (EMT) is a process that the majority of tumors have gone through during tumor progression. The role of EMT in tumorigenesis has been extensively investigated in different cancers including glioma (Lee et al., 2006;Phillips et al., 2006;Hugo et al., 2007;Thiery et al., 2009;Zarkoob et al., 2013). Activated EMT has also been found to be associated with the generation of cancer stem cells (Ye and Weinberg, 2015). KRAS gene polymorphisms are associated with the risk of glioma (Guan et al., 2021). Collectively, GSEA results suggested that LGG patients with high-risk scores tend to develop a faster deterioration than patients with the low-risk scores.
Recently, a growing body of studies has demonstrated how E3 ligases function in the tumor microenvironment (Do et al., 2022;Hosein et al., 2022;Iwamoto et al., 2022). The tumor microenvironment consists of different cells, including tumor cells, tumor stem cells, and stromal cells. These cells form a complex network and interact with one another to regulate tumor malignant behaviors and treatment resistance. Growth factors, chemokines, and cytokines released by immune cells are widely involved in tumor progression and therapy responses (Bindea et al., 2013). Here, we assessed the infiltration level of 28 TIME immune cells in the high-and low-risk groups to explore the roles of identified signature genes. By ssGSEA, our results showed that 25 out of 28 immune cells are more abundant in the high-risk group than in the low-risk group. In line with other studies, macrophages seem to constitute a majority of infiltration in low-grade gliomas (Rossi et al., 1988;Mieczkowski et al., 2015). These could be important findings since it has been shown that the infiltration of macrophages is highly associated with shorter overall survival in low-grade glioma (Müller et al., 2017). Our study observed the increase of most immune cells infiltrated in the high-risk group, which might be correlated with poor prognosis. The immune activity and stromal activity were remarkably elevated in the high-risk group, which once again reiterated the heterogeneity of the glioma TIME between the two groups. Determining the roles of the signature genes in TIME cell infiltration heterogeneity will be beneficial to better understand the mechanisms of the TIME antitumor immune response and developing novel immunotherapy strategies (Kim et al., 2022;Ogino et al., 2022;Tian et al., 2022).
Previous studies have identified different immune checkpoint molecules in gliomas, such as CTLA-4, TIM-3, PD-1, CD48, and LAG3 (Chouaib, 2020). Immunotherapytargeting immune checkpoint proteins have been found to trigger an antitumor immune response (Zeng et al., 2013;Reardon et al., 2016;Boussiotis and Charest, 2018). Therefore, we evaluated the correlation between the risk score and the expression of immune checkpoint molecules. It is worth noting that the risk scores were positively correlated with the expression of immune checkpoint molecules. These results implied the involvement of the signature genes in the pathogenesis of the gliomas via regulating immune-related pathways. Accordingly, a preliminary function analysis for the predicted substrates of the risk signature was used. The results suggested that the substrates might regulate the pathways of cell death and immune responses. This partially explained better responses of individuals with high-risk scores to immunotherapies. Yet, FIGURE 10 | Role of the risk signature in the TIME. (A) ssGSEA scores of the 28 immune cells in the high-risk group and the low-risk group (*p < 0.05; ****p < 0.0001; ns, not significant). (B) Heatmap of the expression of the immune cell markers in the low-risk and the high-risk groups. (C,D) Differences in overall immune and stromal activity between the high-risk group and the low-risk group using the ESTIMATE algorithm. (E) Correlation between the risk signature and 28 TIME cell infiltration. (F) Correlation between the risk signature and immune checkpoints.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 905047 these results further verified the reliability of the risk model in predicting LGG prognosis.
Here, we constructed a prognostic model based on E3RGs and a relevant nomogram in LGG. Of note, we identified a risk signature in the TCGA dataset, which consisted predominately of Caucasian and African American cases. Then, we validated the effect of the risk signature in the CGGA dataset which consisted of Chinese patients. The similar survival analysis results observed in both training and validation sets proved that our model could predict LGG prognosis in varying ethnicities. Therefore, the predictive capability of this model could be beneficial to clinical decision-making with LGG patients. Nevertheless, the construction and validation of the risk model were accomplished by retrospective analysis.
Prospective clinical research needs to be rendered for further validation of this model. Moreover, the molecular mechanism of the genes in the risk model requires in-depth investigation in the future.

CONCLUSION
In summary, by differential expression analyses following univariate Cox regression analysis and log-rank statistical test, E3-related DEGs with a prognostic value were identified. A sevengene risk signature was constructed using the LASSO-Cox regression model. The risk signature achieved good performance in predicting the prognosis of LGG. Patients with FIGURE 11 | PPI network and functional analysis of substrates of the E3RG signature in LGG. (A) PPI network of substrates of the E3RG signature obtained with interaction scores > 0.4 based on the STRING online database and visualized using Cytoscape. (B) Top 10 hub genes identified using the MCC algorithm of the cytoHubba plugin in Cytoscape. (C) GO analysis of substrates of the E3RG signature with three terms biological processes, cellular components, and molecular functions (P.adjust < 0.05). (D) KEGG pathways enriched in the substrates of the E3RG signature analyzed using ClueGO plugin in Cytoscape (p < 0.01).
high-risk scores were more likely to have a poor prognosis compared with the ones with low-risk scores. Functional analyses of the risk signature were performed. This study is expected to benefit the diagnosis and the potential therapeutics of low-grade 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.

AUTHOR CONTRIBUTIONS
ST and PW conceived and designed the study. ST and JZ performed the literature review and conceptualization. ST and PW performed the statistical analyses and wrote the original draft. XS and RS revised the manuscript. All authors read and approved the final manuscript.

FUNDING
This study was supported by grants from the National Natural Science Foundation of China (Grant Number: 81771155). This study was approved by the Ethics Committee of Qilu Hospital of Shandong University.