Abstract
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.
Keywords: E3-related genes, prognosis, tumor immune microenvironment,
risk signature, low-grade gliomas
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
([34]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 ([35]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
([36]Brat et al., 2015; [37]Xia et al., 2018; [38]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 ([39]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 ATP-dependent
reversible process mediated by the ubiquitin-proteasome system (UPS),
including E1 ubiquitin–activating enzymes, E2 ubiquitin–conjugating
enzymes, E3 ubiquitin-protein ligases, and deubiquitinating enzymes
(DUBs) ([40]Reyes-Turcu et al., 2009; [41]Schulman and Harper, 2009;
[42]Buetow and Huang, 2016; [43]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 ([44]Deng et al., 2020). Accordingly,
abnormal ubiquitination may contribute to various human pathologies
such as neurodegeneration ([45]Stieren et al., 2011; [46]Popovic et
al., 2014), autoimmune responses ([47]Zangiabadi and Abdul-Sater,
2022), and oncogenic processes ([48]Rape, 2018). In the UPS, E3
ubiquitin ligase serves as the essential part of the ubiquitination
process owing to its substrate specificity ([49]Zheng and Shabek,
2017). UPS is stringently regulated by E3 ligases that convey the
specificity of ubiquitination. In particular, ubiquitin molecules are
transferred from ubiquitin-conjugating 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 ([50]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 ([51]Que et al., 2022). The degradation of
TUSC2 induced by NEDD4 facilitates glioblastoma progression ([52]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 ([53]Veeriah
et al., 2010; [54]Lin et al., 2015). Tripartite motif-containing
protein 11 (TRIM11), overexpressed in glioma, promotes proliferation,
invasion, migration, and glial tumor growth via the induction of EGFR
([55]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 ([56]Warfel et al., 2011), TRAF6
([57]Feng et al., 2014), and TRIM21 ([58]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)
([59]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 ([60]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.
Materials and Methods
Datasets and E3-Related Gene Acquisition
The transcriptome sequencing data and corresponding clinical data of
primary LGG were obtained from the TCGA database
([61]https://portal.gdc.cancer.gov/) and the Chinese Glioma Genome
Atlas (CGGA) database, respectively, ([62]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
[63]GSE68848 and [64]GSE4290 datasets procured from the Gene Expression
Omnibus database (GEO; [65]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
([66]https://esbl.nhlbi.nih.gov/Databases/KSBP2/Targets/Lists/E3-ligase
s/) and iUUCD2.0 database ([67]http://iuucd.biocuckoo.org/)
([68]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 ([69]Robinson et al.,
2010; [70]Love et al., 2014; [71]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 E3-related
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
([72]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:
[MATH:
Risk score=
∑i=1<
/mn>nCoefi×Expri, :MATH]
(where
[MATH: Coefi
:MATH]
is the coefficient and
[MATH: Expri
:MATH]
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” ([73]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 ([74]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 ([75]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 ([76]Yoshihara et al., 2013). The correlation analysis was
performed based on Pearson’s correlation coefficient and presented
using R packages “corrplot” and “circlize” ([77]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
([78]https://cn.string-db.org/) and visualized using Cytoscape (version
3.9.1) ([79]Shannon et al., 2003; [80]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.
Results
Identification of Differentially Expressed E3-Related Genes With the
Prognostic Value in Low-Grade Gliomas
The overall study workflow is presented in [81]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
([82]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 ( [83]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
([84]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 ([85]Figure 2C).
FIGURE 1.
FIGURE 1
[86]Open in a new tab
Workflow of the study.
FIGURE 2.
[87]FIGURE 2
[88]Open in a new tab
Identification of E3-related DEGs with prognostic value in LGG. (A)
Volcano plot of E3-related DEGs between 451 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.
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
([89]Liu et al., 2021; [90]2022a; [91]2022b; [92]2022c), 20 of which
were selected for further multivariate Cox regression analysis
([93]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
([94]Figure 3C). The risk score was calculated for each LGG patient by
the following formula:
FIGURE 3.
[95]FIGURE 3
[96]Open in a new tab
Construction of the prognostic E3-related signature. (A) LASSO
regression cross-validation for tuning parameter lambda selection. (B)
Coefficient profiles of the LASSO regression analysis. (C) Seven
optimal prognostic E3-related DEGs selected by multivariate Cox
regression analysis are shown in the forest plot. (D) Expression of the
signature genes in the [97]GSE4290 dataset. (E) Expression of the
signature genes in the [98]GSE68848 dataset. (*p < 0.05; **p < 0.01;
***p < 0.001; ****p < 0.0001; ns, not significant).
Riskscore=(0.4121493)∗ AURKA+ (-0.8124184)∗PCGF2+ (0.6318466)∗MAP3K1+
(0.2558092)∗TRIM34+(-0.5951688)∗PRKN+(-0.5669701)∗TLE3+(0.1586579)∗TRIM
17.
The expression of the signature genes was validated in both [99]GSE4290
and [100]GSE68848 datasets and proved consistent with that in the TCGA
cohort ([101]Figures 3D,E). Univariate Cox regression analysis revealed
that seven signature genes were all strongly correlated with the
prognosis ([102]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 ([103]Figure 4B).
FIGURE 4.
[104]FIGURE 4
[105]Open in a new tab
Prognostic analysis of the seven-E3RG risk signature. Univariate Cox
regression analysis (A) and Kaplan–Meier curves (B) of the seven-E3RG
prognostic signature (**p < 0.01; ***p < 0.001).
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 ([106]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 expression of AURKA, MAP3K1, and
TRIM34 was lower in low-risk patients than in high-risk patients, while
PCGF2, PRKN, TLE3, and TRIM17 expressions were higher in patients with
low-risk scores than in those with high-risk scores. The distribution
of risk score, survival status, and the expression of the signature
genes is shown in [107]Figures 5B–D. To evaluate the predictive effect
of the prognostic model, 1-year, 3-year, and 5-year time-dependent ROC
curves were plotted, and the concordance index was calculated. The area
under the curve (AUC) values were 0.9 (1-year ROC), 0.89 (3-year ROC),
and 0.85 (5-year ROC) ([108]Figure 5E). Given the results earlier, the
risk signature presented a superior predictive capacity for LGGs.
FIGURE 5.
[109]FIGURE 5
[110]Open in a new tab
Distribution and prognostic analysis of the E3RG prognostic signature
in the TCGA cohort. (A) Kaplan–Meier curves of overall survival for the
patients in the high-risk group and the low-risk group of the TCGA
cohort. The distribution plots (B), corresponding survival status (C),
and the expression of the seven signature genes (D) in the TCGA cohort.
(E) Time-dependent ROC curve validation of the predictive capacity of
the risk signature in the TCGA cohort.
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 ([111]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 ([112]Figure 6A, p < 0.0001). The distribution of
risk score and survival status is shown in [113]Figures 6B,C. A heatmap
of gene expression in the CGGA cohort is presented in [114]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 ([115]Figure 6E).
Prognostic analyses showed similar results. These results demonstrated
that the E3RG risk signature was positively correlated with LGG
prognosis.
FIGURE 6.
[116]FIGURE 6
[117]Open in a new tab
Validation of the E3RG prognostic signature in the CGGA cohort. (A)
Kaplan–Meier curves of the overall survival in the high-risk group and
low-risk group of the CGGA cohort. The distribution plots (B), survival
status (C), and the expression of seven signature genes (D) in the CGGA
cohort. (E) Time-dependent ROC curve validation of the predictive
capacity of the risk signature in the CGGA cohort.
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 ([118]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
([119]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 independent prognostic factor in
the CGGA cohort ([120]Figures 7C,D). The subtle difference may be
attributed to the ethnicity difference.
FIGURE 7.
[121]FIGURE 7
[122]Open in a new tab
Construction and evaluation of a nomogram. Independent prognostic
factors were identified by univariate and multivariate Cox regression
analyses 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.
On the basis of the seven-E3RG risk signature, we established a
nomogram that could predict the prognosis of LGG in the TCGA cohort
([123]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 ([124]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 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
([125]Figure 8A). As shown by the Venn diagrams in [126]Figure 8B, 528
upregulated genes and 134 downregulated genes in the high-risk group
were identified and applied for GO and KEGG pathway analyses
([127]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
([128]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 ([129]Figure 8D). This
was also consistent with the results of GO analysis.
FIGURE 8.
[130]FIGURE 8
[131]Open in a new tab
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
risk-related 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).
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 ([132]Figure 9).
FIGURE 9.
[133]FIGURE 9
[134]Open in a new tab
Gene set enrichment analysis (GSEA) of DEGs in the high-risk group.
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 TIME-infiltrating
cell types in the low-risk and high-risk groups by ssGSEA
([135]Charoentong et al., 2017). In total, 25 TIME cell types presented
significant differences in infiltration between low-risk and high-risk
groups ([136]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
([137]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 ([138]Figures 10C,D). Correlation analysis was performed, and
potential relations between the risk score signature and each TIME cell
type are shown in [139]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 ([140]Figure 10F), implying the meaningful
roles of risk score signature in predicting the possible response of
LGG patients to clinical immunotherapy.
FIGURE 10.
FIGURE 10
[141]Open in a new tab
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.
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 E3-related gene
signature were searched on Ubibrowser 2.0
([142]http://ubibrowser.bio-it.cn/ubibrowser_v3/). The known substrates
were applied for protein–protein interaction network analysis using the
STRING database ([143]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 ([144]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 ([145]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) ([146]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) ([147]Figure 11D).
FIGURE 11.
[148]FIGURE 11
[149]Open in a new tab
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).
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
([150]Brennan et al., 2013). EGFR could be downregulated by PARK2 but
increased by TRIM11; both are E3 ligases possessing opposite effects on
EGFR ([151]Di et al., 2013; [152]Lin et al., 2015). The PI3K/Akt
pathway is altered in approximately 90% of GBM cases ([153]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 ([154]Winston et al., 1999; [155]Li et al., 2009; [156]Yang et
al., 2009; [157]Chan et al., 2012; [158]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 ([159]Wang et al., 2018;
[160]Nguyen et al., 2021). MAP3K1 might promote glioma stem cell
progression and be positively associated with resistance to
temozolomide (TMZ) and radiotherapy ([161]Bi et al., 2020; [162]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 ([163]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 EGF-stimulated pathways
([164]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 ([165]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 ([166]Wang et al., 2009; [167]Guo et al., 2010;
[168]Zhang et al., 2010). TLE3 expression is positively correlated with
taxane sensitivity in patients with ovarian carcinoma but not breast
cancer ([169]Samimi et al., 2012; [170]Bartlett et al., 2015). TRIM17
augments BRAF-targeted therapy sensitivity of melanoma cells by
preventing BCL2A1 from being ubiquitinated and degraded by TRIM28
([171]Lionnard et al., 2019). The studies mentioned earlier once again
address the potent roles of our risk genes in modulating cancer
biological functions.
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 ([172]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 T cells. Both clinical and
transcriptomic data have uncovered that tumoral MHCII is tightly
correlated with poor prognosis and immune responses ([173]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
([174]Grivennikov et al., 2010). Meanwhile, the inflammatory nature of
the tumor microenvironment promotes the development and survival of
tumors ([175]Mantovani et al., 2008). IL-2/STAT5 signaling is crucial
for the regulation of regulatory T cells and immune tolerance ([176]Shi
et al., 2018). Similarly, the IL6-JAK-STAT3 signaling pathway may
promote tumor cell proliferation, invasion, and metastasis and suppress
immune response ([177]Johnson et al., 2018). In addition, inflammatory
responses and complement-related pathways could affect the
tumorigenesis, immune surveillance, and immunotherapy response
([178]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 ([179]Cirone, 2021;
[180]Pungsrinont et al., 2021). Mutations within genes of this
signaling pathway are the most common events occurring in solid
malignancies including glioma ([181]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 ([182]Unni and Arteaga, 2019). Furthermore,
AKT/PI3K signaling could indirectly activate mTORC1 by the
phosphorylation of PRAS40, a known mTORC1 inhibitor ([183]Sancak et
al., 2007; [184]Thedieck et al., 2007; [185]Vander Haar et al., 2007;
[186]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 ([187]Pusapati et al.,
2016). Recently, treatments targeting PI3K/AKT/mTOR and mTORC1 pathways
have emerged as promising strategies in cancer therapeutics ([188]Yang
et al., 2019; [189]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
([190]Lee et al., 2006; [191]Phillips et al., 2006; [192]Hugo et al.,
2007; [193]Thiery et al., 2009; [194]Zarkoob et al., 2013). Activated
EMT has also been found to be associated with the generation of cancer
stem cells ([195]Ye and Weinberg, 2015). KRAS gene polymorphisms are
associated with the risk of glioma ([196]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 ([197]Do et al., 2022;
[198]Hosein et al., 2022; [199]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 ([200]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 ([201]Rossi et al., 1988; [202]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 ([203]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 ([204]Kim
et al., 2022; [205]Ogino et al., 2022; [206]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 ([207]Chouaib,
2020). Immunotherapy-targeting immune checkpoint proteins have been
found to trigger an antitumor immune response ([208]Zeng et al., 2013;
[209]Reardon et al., 2016; [210]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,
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 seven-gene risk signature
was constructed using the LASSO-Cox regression model. The risk
signature achieved good performance in predicting the prognosis of LGG.
Patients with 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/[211]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.
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.
Publisher’s Note
All claims expressed in this article are solely those of the authors
and do not necessarily represent those of their affiliated
organizations, or those of the publisher, the editors, and the
reviewers. Any product that may be evaluated in this article, or claim
that may be made by its manufacturer, is not guaranteed or endorsed by
the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at:
[212]https://www.frontiersin.org/articles/10.3389/fgene.2022.905047/ful
l#supplementary-material
[213]Click here for additional data file.^ (105.1KB, XLSX)
Abbreviations
WHO, World Health Organization; LGG, low-grade glioma; GBM,
glioblastoma multiforme; PTMs, posttranslational modifications; UPS,
ubiquitin-proteasome system; DUBs, deubiquitinating enzymes; E3RGs,
E3-related genes; PROTAC, proteolysis targeting chimeric; TCGA, The
Cancer Genome Atlas; CPM, counts per million reads mapped; FPKM,
Fragments Per Kilobase of transcript per Million mapped reads; DEGs,
differentially expressed genes; LASSO, least absolute shrinkage and
selection operator; HRs, hazard ratios; CIs, confidence intervals; ROC,
receiver operating characteristic; GO, Gene Ontology; KEGG, Kyoto
Encyclopedia of Genes and Genomes; BP, biological process; CC, cellular
component; MF, molecular function; GSEA, gene set enrichment analysis;
ssGSEA, single-sample gene set enrichment analysis; PPI,
protein–protein interaction network; AUC, area under the curve; TIME,
tumor immune microenvironment; GSCs, glioma stem cells; TMZ,
temozolomide; MHC, major histocompatibility complex; ECM, extracellular
matrix.
References