Function of Differentially Expressed Genes and Risk Genes in Glioma
Identifying the relatively typical biological processes in gliomas is the basis for seeking markers. Therefore, we first conducted differential gene enrichment analysis (KEGG and GO) based on glioma (containing the glioblastoma multiforme and low-grade glioma) sample data from TCGA database. Based on significantly upregulated genes, the enrichment results of both KEGG and GO reflected entries related to cell cycle regulation. We combined with the characteristics of tumor cell cycle, it was suggested that this process is the basis of malignant phenotype or progression of tumor cells. (Fig. 1A, B, where the red arrow is pointing). Meanwhile, there were no typical characteristic items in the significantly down-regulated gene enrichment results (Supplementary Fig. 1A, B). To further explore the commonality between the function of risk genes and up-regulated genes, 8344 risk genes (HR > 1, P < 0.05) were obtained by analyzing TCGA-glioma data, and they were compared with 1091 up-regulated genes (Log2FC > 1, P < 0.05) for Venn diagram analysis. The results showed that there were 679 genes with common characteristics (Fig. 1C). Enrichment analysis of these genes revealed four biological processes associated with the cell cycle in the top 5 entries (Fig. 1D, where the red arrow is pointing). These results suggest that upregulated risk genes in gliomas may also be involved in cell cycle processes.
Based on the above enrichment results, we conducted Venn diagram analysis for the genes from these four items, which found that 27 genes were the common part (Fig. 1E). Correlation analysis of these genes in glioma showed that 17 genes were strongly correlated with each other (R > 0.8, P < 0.05), and all showed high risk at high expression levels (Fig. 1F, G). These results confirmed the high-risk characteristics for the selected genes. Meanwhile, Cox regression analysis showed that 17 genes were significant risk factors in univariate analysis, and the multivariate prognostic analysis shows that CHEK1, UBE2C and BUB1 were independent prognostic risk factors (Supplementary Fig. 1C, D, HR > 1, P < 0.05). In addition, according to the expression level of 17 genes in glioma, UBE2C had the most significant upregulation (Fig. 1H, ***P < 0.001). In summary, UBE2C may be a risk factor and significantly upregulated in gliomas, which may be the internal cause of malignant phenotype of tumors.
Clinical Classification and Prognostic Significance of UBE2C in Glioma
UBE2C was selected as a typical representative gene. To further verify its clinical significance in glioma, these samples from the TCGA database were clinically classified (Supplementary Table 1). The results showed that the expression level of UBE2C was significantly up regulated with the progression of the disease in both WHO classification (G2, G3, G4) and primary therapy outcome (CR, Complete response; PR, Partial Response; SD, Stable disease; PD, Progressive disease) (Fig. 2A, B, **P < 0.01, ***P < 0.001). Among the pathological grades of glioma, the malignant degree of glioblastoma was the highest [5, 19], and the expression level of UBE2C was significantly higher than that of the other three types (Fig. 2C, ***P < 0.001). In terms of age and survival status, glioma patients over 60 years of age and those who had died showed significant UBE2C upregulation (Fig. 2D, E, ***P < 0.001). These features can also be demonstrated by the results of Cox regression analysis of clinical grading (Supplementary Table 2). These results suggest that the expression level of UBE2C is significantly different in the clinical grades of glioma, which can significantly affect the prognosis of these clinical grades.
In addition, the ROC Curve (Receiver Operating Characteristic Curve) analysis of UBE2C in glioma samples showed that the AUC (Area Under ROC Curve, CI: 0.935–0.993) value was 0.964. This result indicates that UBE2C gene has good diagnostic value for glioma (Fig. 2F). To assess the prognostic value of UBE2C in glioma, survival curves (K-M) with OS (Overall Survival), DSS (Disease Specific Survival), and PFI (Progress Free Interval) were performed on these clinical samples, respectively. The results showed that the three types of K-M curves showed that the survival rate of the group with high expression of UBE2C was relatively low, that is, the expression level of this gene was significantly negatively correlated with the survival rate of patients (Fig. 2G-I). In conclusion, the differences in the expression of UBE2C in clinical classification and survival prognosis can reflect the related characteristics of disease progression, and thus show certain clinical diagnosis and prognostic value.
The Intrinsic Effect of UBE2C in Glioma
Based on the clinical significance of UBE2C in gliomas, we need to further explore the internal influence of this gene in gliomas, which can reflect the value of UBE2C as a marker of endogenous adverse factors. We need to further explore the major biological processes involved in UBE2C by enrichment analyzing the genes that were significantly positively correlated with UBE2C (R > 0.6, adj. P < 0.05). The results show that the TOP 20 items presented mainly include two biological processes, which are cell cycle regulation (red arrow) and DNA behavior (blue arrow) (Fig. 3A). The results of corresponding network and correlation analysis also reflect a fact that the cell cycle and DNA behavior are two major units, which is a close bond between them (Fig. 3B, C, red and blue dashed line ellipses). GSEA analysis of differential genes based on TCGA verified the significance of cell cycle and DNA behavior, respectively (NES > 1, adj. P < 0.05, FDR q < 0.25) (Fig. 3D, E). The analysis results of the GEO-sourced dataset (GSE12657) can also fully verify this conclusion (Supplementary Fig. 2A, B, NES > 1, adj. P < 0.05, FDR q < 0.25). Therefore, the above results indicate that UBE2C may satisfy the adaptive survival of tumor cells in glioma mainly via the regulation of cell cycle and DNA behavior.
Summarizing the above results, it can be inferred that UBE2C gene may be related to malignant proliferation of glioma cells. To verify this association, we separately display the IHC (Immunohistochemistry) images from The Human Protein Atlas (THPA) database. Since Ki67, PCNA and MCM7 are common markers reflecting tumor cell proliferation [20], they are related to the replication behavior of nuclear DNA. Therefore, according to the IHC of these three proliferative markers, the nuclear coloring degree of tumor region is deeper than adjacent region (Fig. 4B-D). Meanwhile, the staining trend of UBE2C cells in tumor region was consistent with these three proliferative markers (Fig. 4A). In addition, the scatter plot of correlation analysis shows that UBE2C had a significant positive correlation with Ki67, PCNA and MCM7, respectively (Fig. 4E-G). In summary, the expression level of UBE2C is related to the proliferation of tumor cells, and it may be an internal factor affecting the adaptive survival of tumor cells.
To verify the expression of UBE2C in glioma cell lines, our analysis based on the mRNA data collected in the CCLE database (Cancer Cell Line Encyclopedia), which showed that there was little difference in the expression of UBE2C in 48 brain cancer cell lines (Fig. 4H). Meanwhile, the differences between the common glioma cell lines LN229 and U251MG were relatively obvious (Fig. 4H, where the red arrow is pointing). The expression of UBE2C in 4 glioma cell lines (LN229, U87MG, A172 and U251MG) was detected by qPCR and Western blot, respectively. The results showed that the mRNA expression level of UBEC in these four cell lines was significantly higher than that of normal glioma cells (HEB), but there was no difference among the four cells (Fig. 4I, **P < 0.01, ***P < 0.001). In addition, the expression of UBE2C protein in the four types of glioma cells was similar to mRNA levels (Fig. 4J, *P < 0.05, **P < 0.01, ***P < 0.001). In conclusion, the relatively high expression level of UBE2C in glioma cell lines was consistent with TCGA data and IHC results.
Effect of UBE2C Expression on Malignant Phenotype for Glioma Cell
Based on the difference of UBE2C expression in glioma (tissue and cellular level) and its relevance to the clinical prognosis of this tumor, we need to further verify the effects of UBE2C on the malignant phenotype for the glioma at the cellular level. UBE2C was highly expressed in glioma cell lines LN229 and U251MG (Fig. 4H, J). Therefore, we achieved knockdown of UBE2C protein expression via transfecting these two cells with siRNA. The results showed that siRNA3 had the most significant knockdown effect on these two types of cells (Fig. 5A, B, *P < 0.05, **P < 0.01). Since the knockdown effects of siRNA3 above, cell scratch experiments showed that the knockdown of UBE2C in both LN229 and U251MG cells could significantly inhibit the healing ability of these two cells (Fig. 5C, D, **P < 0.01, ***P < 0.001). In addition, when UBE2C is knocked down in LN229 and U251MG cells, it can significantly inhibit the migration ability of these two cells (Fig. 5E, F, **P < 0.01). These results indicate that the expression level of UBE2C is positively correlated with the migration ability of tumor cell lines LN229 and U251MG.
In addition, the cell cloning assay of LN229 and U251MG showed that when UBE2C was knocked down, it could significantly inhibit the ability of clonal formation for the two tumor cells (Fig. 5G, H, **P < 0.01, ***P < 0.001). Meanwhile, the cell proliferation curve assay (OD450, CCK8) showed that the knockdown of UBE2C could significantly inhibit the proliferation ability of these two types of cells (Fig. 5I, J, *P < 0.05, **P < 0.01). Therefore, these results indicate that the expression level of UBE2C can significantly affect the malignant phenotype of LN229 and U251MG cells.
The External Effects of UBE2C in Glioma
As an intrinsic risk factor affecting the progression of glioma, UBE2C may be involved in the regulation of cell cycle to positively link the proliferative activity of tumor cells. However, we also need to understand its correlation with external factors. Based on TCGA data, we analyzed the correlation between UBE2C and 24 types of immune cell infiltration in glioma. The results showed that the expression level of UBE2C was most significantly correlated with the immune infiltration of Th2 cells (Fig. 6A, R = 0.871, ***P < 0.001). Therefore, we are concerned that this feature may be the focus as UBE2C is linked to external factors of the tumor. According to the secretory phenotype of Th2 cells [21, 22] and the trend of UBE2C expression, and combined with their correlation analysis, we can see that IL6 (R = 0.396, P < 0.01) and IL10 (R = 0.338, P < 0.01) are the most significant (Fig. 6B, C). In conclusion, UBE2C may be closely related to the immune infiltration of Th2 cells.
It has been reported that IL6 or IL10 can induce immunosuppressive phenotypes in the tumor microenvironment, which leads to the immune escape mechanism of tumor cells [23]. Therefore, for malignant progression and poor prognosis of glioma, the correlation between UBE2C and external adverse factors may reflect the essence behind this phenomenon. We analyzed the correlation between the expression of UBE2C and 14 immune checkpoints in glioma, and the results showed that eight immune checkpoints (CD80, SIGLEC7, LAG3, CD28, PDCD1LG2, PDCD1, HAVCR2, SIGLEC15) were significantly positively correlated with UBE2C (Fig. 6D, E, R > 0.3, P < 0.001). Based on the correlation analysis of Th2 cells secretion phenotype and immune checkpoints, we found that IL6 and IL10 were most significantly correlated with various immune checkpoints, among which six immune checkpoints were more typical (where the red arrow is pointing) (Fig. 6F, R > 0.5, *P < 0.05). These results suggest that secretory phenotypes IL6 and IL10 from Th2 cells may be important factors affecting immunosuppressive phenotypes in gliomas.
Based on the systematic correlation characteristics of UBE2C - IL6/IL10 - immune checkpoint axis, we further analyzed the prognostic value of glioma. The results showed that the expression levels of 13 genes, such as IL6/IL10 and immune checkpoint, were significantly negatively correlated with the survival rate of patients (Fig. 7A). Meanwhile, IHC images from the THPA database displayed relatively typical expressions of six proteins, which indicated that IL6, IL10, PDCD1LG2, PDCD1, HAVCR2 and CD80 were upregulated in glioma tissues, respectively (Fig. 7B-G). Taking together, the expression levels of UBE2C-related immune infiltrating Th2 cells and immune checkpoints can significantly affect the survival and prognosis of glioma, which suggesting that UBE2C may be an external risk factor to predict tumor progression.
Correction of UBE2C-Related Risk Factors with Cell Invasion in Glioma
Since the high aggressiveness of glioma cells has been a major difficulty in the treatment of this disease [24, 25], we should also consider its association with cell aggressiveness phenotype when exploring the internal and external correlation of UBE2C. It has been reported that the invasiveness of tumor cells is correlated with soluble immune checkpoints, which reflects that immune checkpoints in the immune microenvironment can significantly affects the invasive phenotype of tumor cells [26]. Therefore, we performed network and correlation analyses of between MMPs (Matrix Metalloproteinase) family genes and immune checkpoints, which due to the MMP family plays an important role in the invasion of tumor cells [27]. Based on these backgrounds and analysis, we found that UBE2C-related risk genes (IL6, IL10 and immune checkpoints) were significantly positively correlated with MMP1, -2, -3, -7, -8, -9, -10, -11, -12, -13, -14 and − 19, respectively (Fig. 8A, *P < 0.05). In addition, online analysis by STRING showed that UBE2C, IL6, IL10 and MMP9 acted as a “linker” that between immune checkpoints and the MMP family (Fig. 8B). Subsequently, after verifying the correlation between seven MMPs (The most typical seven MMPs, which including MMP1, -2, -7, -9, -11, -14, -19) and immune infiltration in glioma, it was found that MMP2, MMP9 and MMP11 had the most significant correlation with Th2 cells in glioma (R > 0.5, P < 0.001) (Fig. 8C-E, Supplementary Fig. 3A-D). In conclusion, UBE2C, IL6 and IL10 may be the link between immune checkpoint and MMP family. This result also confirmed that the immune invasive phenotype (IL6, IL10 and immune checkpoints) of Th2 may be related to the aggressiveness of tumor cells.
Based on these findings, we need to further clarify the prognosis of these seven MMPs for the glioma patients. K-M survival analysis (OS) showed that the expression levels of seven MMPs were significantly negatively correlated with the survival rate of glioma patients, and the MMP2, MMP9 and MMP11 were more typical (Fig. 8F-H, Supplementary Fig. 3E-H). IHC images from the THPA database showed that MMP2 and MMP14 expression differences were most significant in gliomas (MMP1 and MMP19 data were not included) (Supplementary Fig. 4A-E). Finally, multivariate Cox regression models of UBE2C, IL6, IL10, MMP2, MMP9 and MMP11 were analyzed by nomogram to predict 1-, 3-, and 5-year survival probabilities. The results showed that the contribution value of UBE2C, IL10 and MMP2 was the largest (total score ratio: 250/280 = 89.3%), which may be used as a more accurate multi-factor prediction model for glioma patients, and the prediction model was consistent with the conclusion, that is, UBE2C was correlated with immune invasion and cell invasion (Fig. 8I). Meanwhile, the calibration curve verifies that the 1-, 3-, and 5-year survival probability curves fit well with the ideal line based on the same conditions, which indicates that the multi-factor prediction model is reliable (Fig. 8J). In conclusion, among the risk factors for glioma, UBE2C may be systematically associated with IL10 and MMP2, which become a better predictive model. Therefore, UBE2C may be an ideal biomarker for the prognosis of glial patients.