Construction of Novel Methylation-Driven Gene Model and Investigation of PARVB Function in Glioblastoma

Background Glioblastoma multiforme (GBM) is characterized by widespread genetic and transcriptional heterogeneity. Aberrant DNA methylation plays a vital role in GBM progression by regulating gene expression. However, little is known about the role of methylation and its association with prognosis in GBM. Our aim was to explore DNA methylation-driven genes (DMDGs) and provide evidence for survival prediction and individualized treatment of GBM patients. Methods Use of the MethylMix R package identified DMDGs in GBM. The prognostic signature of DMDGs based on the risk score was constructed by multivariate Cox regression analysis. Receiver operating characteristics (ROC) curve and C-index were applied to assess the predictive performance of the DMDG prognostic signature. The predictive ability of the multigene signature model was validated in TCGA and CGGA cohorts. Finally, the role of DMDG β-Parvin (PARVB) was explored in vitro. Results The prognostic signature of DMDGs was constructed based on six genes (MDK, NMNAT3, PDPN, PARVB, SERPINB1, and UPP1). The low-risk cohort had significantly better survival than the high-risk cohort (p < 0.001). The area under the curve of the ROC of the six-gene signature was 0.832, 0.927, and 0.980 within 1, 2, and 3 years, respectively. The C-index of 0.704 indicated superior specificity and sensitivity. The six-gene model has been demonstrated to be an independent prognostic factor for GBM. In addition, joint survival analysis indicated that the MDK, NMNAT3, PARVB, SERPINB1, and UPP1 genes were significantly associated with prognosis and therapeutic targets for GBM. Importantly, our DMDG prognostic model was more suitable and accurate for low-grade gliomas. Finally, we verified that PARVB induced epithelial-mesenchymal transition partially through the JAK2/STAT3 pathway, which in turn promoted GBM cell proliferation, migration, and invasion. Conclusion This study demonstrated the potential value of the prognostic signature of DMDGs and provided important bioinformatic and potential therapeutic target data to facilitate individualized treatment for GBM, and to elucidate the specific mechanism by which PARVB promotes GBM progression.


INTRODUCTION
Glioblastoma multiforme (GBM) is the most frequent primary brain tumor in adults. Despite aggressive treatment regimens, GBM patients generally have a poor prognosis of less than 14 months (1). Rapid development in bioinformatics has greatly contributed to explorations of the molecular characteristics of cancer. As well, specific GBM molecular markers have been discovered, providing new insights into the progression mechanism, diagnosis, and treatment of GBM (2). Among the many molecular markers, isocitrate dehydrogenase (IDH)1/2 mutation and O 6 -methylguanine DNA methyltransferase (MGMT) methylation are the best-known. Compared with GBM patients without IDH mutations, GBM patients with IDH mutations have longer survival times (3,4). GBM patients with MGMT methylation are more sensitive to temozolomide (TMZ) (5,6). Therefore, it is clinically important to identify effective and promising biomarkers for the prognosis of GBM.
DNA methylation modification is an important part of epigenetics, which contributes to the transcriptional regulation of genes and maintenance of genomic stability (7). DNA methylation has low variability and relative semi-stability, but is closely related to cell process activity. The DNA methylation status of the CpG island of the promoter can regulate the expression of tumor-related genes and plays a key role in the occurrence and development of cancer (8). Moreover, the high plasticity of DNA methylation allows tumor cells to rapidly adapt to metabolic constraints or physiological changes during tumorigenesis (9). Therefore, the combination of transcriptome and methylation status can help identify novel markers, improve cancer diagnosis, and predict clinical outcomes.
b-Parvin (PARVB) is a member of the parvin protein family. The protein localizes to focal adhesions and plays important roles in cell adhesion, proliferation, and migration by activating multiple signaling pathways (10). Overexpression of PARVB has been associated with poor prognosis in human colorectal cancer (11) and tongue squamous cell carcinoma (12). PARVB has also been associated with epithelial-to-mesenchymal transition (EMT), a biological process in which epithelial cells undergo multiple biochemical changes, ultimately switching to a mesenchymal phenotype. Cells undergoing EMT lose their apical basolateral polarization and acquire a fibroblast-like morphology, characterized by weak cell adhesion and enhanced ability of the cell to migrate and spread (13,14). The Janus Kinase/Signal Transducer and Activator of Transcription 3 (JAK/STAT3) signaling pathway is important in various types of cancer. Activation of this pathway leads to increased tumorigenic and metastatic capacity by enhancing EMT (15).
In this study, based on the Wilcoxon rank test, genes with significantly different expression in gliomas were selected from the transcriptome information in The Cancer Genome Atlas (TCGA) and Genotype-Tissue Expression (GTEx) databases. An integrative approach was used to identify GBM-related methylation-driven genes by combining the transcriptome and DNA methylation profiles from the TCGA database. A Cox survival predictive model-based risk score of six DNA-drive methylation genes (DMDGs) was successfully constructed. The score was effective in determining GBM patients with poor prognosis and displayed stronger predictive power in lowgrade glioma (LGG) patients. The verification that EMT induced by PARVB can enhance the proliferation, migration, and invasion of GBM cells through the JAK2/STAT3 pathway emphasizes the value of the PARVB gene as a potential therapeutic target. Our findings provide new insights into the molecular mechanisms of GBM and prompt a more individualized therapy for this prevalent disease.

Selection of Differentially Expressed Genes in Glioma
RNA-seq transcriptome data obtained from normal tissues in GTEx and gliomas in TCGA were standardized using the limma package. The cutoff criteria were set as | log2 fold change (FC) | > 0.5 and p < 0.05. DEGs between normal and glioma tissues were selected based on the Wilcoxon rank test.

Selection of DMDGs
The cutoff criteria of the DEGs and aberrantly methylated genes (AMGs) were set as | log2 FC | > 0.2 and p < 0.05 between LGG (selected astrocytoma only) and GBM samples. Combining the data of DEGs and AMGs in GBM, the MethylMix R package was used to screen DMDGs (16). The correlation coefficient was set to Cor < −0.4. DMDGs meeting the cutoff criteria were screened. The expression and methylation of these DMDGs were visualized based on the heatmap R package.

Functional and Pathway Enrichment Analyses of DMDGs
Gene ontology (GO) analysis of DMDGs was performed using the Database of Annotation, Visualization and Integrated Discovery (DAVID) (https://david.abcc.ncifcrf.gov/). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using KOBAS 3.0 (http:// kobas.cbi.pku.edu.cn/). The GOplot R package was used to visualize the significantly enriched GO terms and KEGG pathways. Additionally, KEGG pathway analysis of PARVB was performed using Gene Set Enrichment Analysis (GSEA) software (www.gsea-msigdb.org). The cutoff criterion was set at p < 0.05.

Construction, Evaluation, and Validation of DMDG Prognostic Signature
Univariable Cox regression analysis was utilized to select prognosis-related DMDGs (threshold value of genes was set as p < 0.05) in GBM patients. Applying multivariate Cox regression analysis optimized the constructed model by calculating the optimal AIC value. The smaller the AIC, the better the model, and the one with the smallest AIC is usually chosen. AIC=(2k-2L)/n. Genes with high correlation were removed in this process, and finally genes with large differences were selected to construct the model. Subsequently, the DMDG prognostic signature was constructed by the linear combination of the expression levels of DMDGs using the b coefficient calculated from multivariate Cox regression as the weight. The risk score for each patient was calculated as A (expression level of gene) *(beta1) + B*(beta2) + C*(beta3) +……. + (N-2)*(beta(n-2))+ (N-1)*(beta(n-1))+ (N)* (beta(n)). By setting the median value of the risk score as the cutoff value, GBM patients were separated into high-risk and low-risk groups. The difference in overall survival (OS) between the two groups was evaluated by Kaplan-Meier survival analysis using the log-rank test. Receiver operating characteristic (ROC) curves based on the survivalROC package were applied to assess the predictive performance of the prognostic signature of DMDGs within 1, 2, and 3 years. The predictive accuracy of the DMDG prognostic signature was evaluated using Harrell's C-index with the survcomp package. The constructed risk model was re-scored by four previously published methylation-based gene signatures. The signature derived in this study was compared with these four other signature ROC curves. In addition to univariate Cox regression analysis, multivariate Cox regression analysis was performed to evaluate the resolution of the DMDG prognostic signature. Finally, the DMDG prognostic signature constructed using the TCGA database was validated in the CGGA database.

Combined Survival Analysis Based on Expression and Methylation of DMDGs
Using the median values of gene expression and DNA methylation levels of DMDGs, patients were divided into low expression patients with hypermethylation and high expression patients with hypomethylation. The survival differences between the two groups were evaluated using Kaplan-Meier survival analysis.

Cell Lines and Cell Culture
The LN229 human GBM cell line and HG7 patient-derived GBM cells were cultured as previously described (17). Cell line authentication was performed using short tandem repeat analysis (Genetic Testing Biotechnology, Jiangsu, China). Cell lines were actively passaged for up to 6 months. Only cells below passage 15 were used for the experiments.

Chemical Reagents, Antibodies, and Transfection
The JAK2/STAT3 inhibitor WP1066 was purchased from Selleck Chemicals (Boston, MA, USA). PARVB overexpression plasmids and short hairpin RNA (shRNA) were purchased from GeneChem (Shanghai, China). Ara-C (cytarabine) was purchased from Sigma-Aldrich (St Louis, MO). Cells were transiently transfected with overexpression plasmids or shRNA for 2 weeks of puromycin selection to obtain stable cells. Fortyeight hours after transfection, western blot analysis of the collected cells was performed.

Immunoblotting and Immunohistochemistry
GBM samples or cells were harvested using RIPA supplied with proteinase and phosphatase inhibitor cocktails (Selleck Chemicals, Shanghai, China). The protein concentration was determined using a BCA Assay Kit. Western blot assays were performed as described previously (18). The antibodies used are listed in Supplementary Table S1.
Paraffin-embedded GBM tissues were obtained from patients treated at First Affiliated Hospital of Nanchang University who provided informed consent. The study was approved by the hospital's institutional ethics committee. Immunohistochemistry was performed as described previously (19). Clinical GBM samples were immunostained with primary antibodies against PARVB, E-cadherin, Claudin-1, N-cadherin, and vimentin at 4°C overnight. After washing in PBS, the tissues were exposed to secondary antibody (horseradish peroxidase [HRP]-conjugated anti-mouse/rabbit polymeric antibody) for 30 min. The 3,3'diaminobenzidine (DAB) Staining Kit (Zsgb Bio Inc., Beijing, China) was used as a chromogen for 1 min of incubation to allow for proper brown color development.

Cell Proliferation Assay
Cell Counting Kit 8 (Dojindo, Kumamoto, Japan) was used to detect cell proliferation. The fluorescence at 450 nm was recorded using a microplate reader.

Colony Formation Assay
After pancreatin digestion of logarithmic growth phase cells, the cells were harvested and resuspended in complete medium [basal medium + 10% fetal bovine serum (FBS)]and counted. Cells were seeded at 800 cells/well in 6-well plate culture plates for each experimental group. Cells were cultured for 8 days. One milliliter of Crystal violet dye solution was added to each well and the cells were stained for 30 min. Cells were washed several times with PBS, the number of colonies was captured using an Olympus camera (Tokyo, Japan) and counted using ImageJ software.

Three-Dimensional (3D) Spheroid Cell Invasion Assay and Transwell Migration Assay
The 3D spheroid cell invasion assay was performed using a Cultrex ® 96-well 3D Spheroid BME Cell Invasion Assay kit (Trevigen, Gaithersburg, MD, USA). Approximately 3×10 3 cells were resuspended in 50 µL of 1× spheroid formation ECM and were added to each well of a 96-well plate (Helgerman CT, Gaithersburg, MD, 20877). The plate was centrifuged at 200 ×g for 5 min, and incubated at 37°C for 72 h to promote the formation of cell spheres. Next, 50 mL of the invasion matrix/well plate was placed on ice for 10 min, followed by centrifugation at 300×g for 5 min at 4°C. The plates were then incubated at 37°C for 3 days. The cells were imaged using an Axio Imager M2 microscope (Carl Zeiss, Jena, Germany).
Transwell assays were performed in 24-well cell culture chambers with 8 mm pore Transwell inserts precoated with Matrigel. In brief, cells were seeded in 200 mL of culture medium containing 1% FBS. Five hundred microliters of medium containing 50% FBS was added to each lower chamber. After 24 h, the cells were fixed and stained with hematoxylin and eosin. Cells passing through the membrane were imaged (10×) and randomly counted into five independent fields.

Statistical Analyses
Statistical analyses were conducted using the R statistical package (R version 4.0.2). A two-tailed p < 0.05 was considered a significant difference between groups. Significant differences between the groups were estimated using the Student's t-test. One-way analysis of variance (one-way ANOVA) was performed for at least three groups.

Selection of Significant DEGs in Glioma
The cutoff criteria were set as | log2 FC | > 0.5 and p < 0.05. Based on RNA-seq data of 1152 normal tissues in GTEx and 697 gliomas (529 LGG and 168 GBM) in TCGA, 7006 DEGs were selected based on the Wilcoxon rank test between normal and glioma tissues. Of these, 3137 genes were upregulated and 3845 genes were downregulated (Supplementary Table S2). A filtering flow chart for the DMDGs is shown in Figure 1.

Identification of DMDGs in GBM
The 7006 DEGs selected between normal and glioma tissues (Supplementary Table S2) were analyzed by the Wilcoxon rank test to further screen for DEGs and AMGs in 195 LGG (astrocytoma only) and 63 GBM (complete transcriptome and DNA methylation profiles) samples in TCGA. In the MethylMix R package, the screening criteria of DEGs and AMGs was set as | log2FC | > 0.2, p < 0.05 and Cor < −0.4. Seventy-two DMDGs were identified, of which nine were hypermethylated and 63 were hypomethylated (Supplementary Table S3). The expression level and methylation value of DMDGs are shown in a heat map (Figure 2A, B and Supplementary Figure S1A).

Functional Enrichment Analysis of DMDGs in GBM
The possible functions and pathways of the 72 DMDGs in GBM were explored using GO functional enrichment analysis and KEGG pathway enrichment analysis. Among biological processes (BP), molecular function (MF), and cellular component (CC), functional analysis showed that DMDGs were enriched in positive regulation of EMT, regulation of cell morphogenesis, mesenchyme development, mesenchymal cell differentiation, positive regulation of inflammatory response, and positive regulation of oligodendrocyte differentiation ( Figure 3A). In addition, KEGG pathway enrichment analysis showed that the DMDGs were enriched in programmed cell death ligand 1 expression and programmed cell death protein 1 (PD-1) checkpoint pathway in cancer, cellular senescence, and advanced glycation end product (AGE)-receptor for AGE (RAGE) signaling pathway in diabetic complication phenotype-related pathways ( Figure 3B).

Evaluation and Comparison of DMDG Prognostic Model
According to the median DMDG risk score, GBM patients with complete clinical data were divided into high-risk and low-risk cohorts. The distribution of risk scores for the cohorts is shown in Figure 6A. As the risk score increased, the prevalence of death increased and the number of surviving patients decreased. The methylation patterns of the six DMDGs in GBM are  shown in Figure 6B. Survival analysis showed that the low-risk cohort had longer OS than the high-risk cohort ( Figure 6C). ROC curve analysis showed that the area under the curve (AUC) values of the ROC curve were 0.832, 0.927, and 0.980 within 1, 2, and 3 years, respectively ( Figure 6D).  Figure 6E). We compared our signature with four previously published GBM methylation-based gene signature panels: the 9gene methylation signature (20), 16-gene methylation signature

Validation of DMDG Prognostic Model
The predictive ability of the 6-gene model was validated in a GBM patient TCGA cohort (n = 88) who were not involved in the construction of the model. Based on the cutoff value (-3.22), 56 patients were assigned to the high-risk group and 32 patients to the low-risk group. The expression pattern of these six DMDG prognostic signatures in the low-and high-risk groups are shown in Figure 7A. The distributions of the risk scores are shown in Figure 7B. Moreover, the distributions of risk scores and OS status of each patient in Figure 7C indicated good discrimination between the low-risk and high-risk groups. Survival curves demonstrated that high-risk patients had a significantly poorer prognosis than low-risk patients (p < 0.01) ( Figure 7D). The AUCs of the ROC curves were 0.607, 0.760, and 0.655 within 1, 2, and 3 years, respectively ( Figure 7E). We also validated the prognostic value of the risk signature in an LGG patient TCGA cohort (n = 529). Based on the cutoff value (-3.22), LGG patients with complete clinical data were divided into high-risk (n = 71) and low-risk cohorts (n = 458). The expression pattern of the prognostic signature, risk scores, OS status, and survival analysis (p < 0.001) are shown in Supplementary Figures S2A-D). The AUCs of the ROC curves were 0.830, 0.774, 0.748, and 0.680 within 1, 2, 3, and 5 years, respectively (Supplementary Figure S2E). In addition, the prognostic power of the risk model was validated in an LGG (astrocytoma only) patient TCGA cohort (n = 195). Based on the same cutoff value, a total of 37 patients were categorized into the high-risk group, while the remaining 158 patients were classified into the low-risk group. The expression pattern of DMDG prognostic signature, risk scores, and OS status are shown in Figure 7F-H. The low-risk column had a more optimized OS than the high-risk column (p < 0.001) ( Figure 7I). The AUCs of the ROC curves were 0.800, 0.834, 0.873, and 0.701 within 1, 2, 3, and 5 years, respectively, demonstrating the excellent predictive ability of the risk score ( Figure 7J). In addition to TCGA database, glioma patients (n = 127) available in the CGGA database were also used to validate our risk model ( Figure 7K-O). The main feature of the present study is that grouping occurred between LGG and GBM, but not between normal and GBM. The study limited the singularity of tumor types, in which all tumors chose astrocytoma (selection astrocytoma only in LGG, while GBM was a high-grade astrocytoma). Therefore, our risk model was more suitable for LGG (especially astrocytoma) and had a higher prediction accuracy.

Overexpression of PARVB Is Correlated With Shorter Survival and Upregulated PARVB-Induced EMT in GBM Patients
We performed a single gene survival analysis for the prognostic model, dividing the genes into low methylation high expression group and high methylation low expression group, using the median of methylation value as the cutoff values. MDK, SERPINB1, NMNAT3, PARVB, and UPP1 had prognostic value. The low methylation, high expression group was associated with poor prognosis (p < 0.05; Supplementary Figure S3A-E). PARVB has a broad spectrum of biological actions during cancer development (11,12), but its mechanistic exploration in GBM has been less studied and was therefore chosen for further study. Immunohistochemical analysis revealed that PARVB expression was positively correlated with the mesenchymal markers N-cadherin and Vimentin, and negatively correlated with the epithelial markers E-cadherin and Claudin-1 in GBM tissues ( Figure 8A). Survival curves demonstrated that high expression PARVB patients had a significantly poorer prognosis than low expression PARVB patients in 20 GBM patients (p = 0.0088) ( Figure 8B).
Western blot analysis showed that PARVB overexpression led to the downregulation of epithelial markers and upregulation of mesenchymal markers, whereas PARVB knockdown resulted in the opposite effect ( Figure 8D). CCK8, colony formation, Transwell, and 3D invasion assays showed that PARVB overexpression promoted proliferation, migration, and invasion  in LN229 and HG7 cells. In contrast, PARVB knockdown significantly inhibited proliferation, migration, and invasion of LN229 and HG7 cells (Figures 8C, E-G). In addition, we observed same effects in the results of migration and invasion even though treated with Ara-C (cytarabine, 3mg/mL), which was usually used to inhibit cell proliferation. (Supplementary  Figures S3F, G). These data suggest that PARVB promotes malignant behavior by regulating EMT in GBM.

PARVB Promotes EMT by Activating the JAK-STAT Pathway
To further investigate the mechanism underlying the promotion of EMT by PARVB, GSEA of single genes was performed. The analysis demonstrated that the high PARVB expression group was significantly enriched in the JAK-STAT signaling pathway (p < 0.001; Figure 9A). Knockdown of PARVB decreased phosphorylated (p)-JAK2 and p-STAT3 levels, and PARVB overexpression increased p-JAK2 and p-STAT3 levels ( Figure 9B). JAK2 inhibitor prevented PARVB-induced EMT in LN229 and HG7 cells ( Figure 9C). To clarify whether the promotion of EMT by PARVB involved the JAK2/STAT3 pathway, the biological effects of the JAK2/STAT3 pathway were assessed in PARVB-overexpressing GBM cells. Inhibition of JAK2/STAT3 signaling prevented PARVB-induced proliferation, migration, and invasion phenotypes in GBM cells overexpressing PARVB ( Figure 9D-F). We also observed similar results of migration and invasion even though treated with Ara-C (Supplementary Figure S4A). These data demonstrate that PARVB promotes EMT through the JAK2/STAT3 pathway.

DISCUSSION
Glioma is the most common malignant brain tumor, and its prognosis is not ideal. GBM is highly malignant with high mortality (24). Even with standard treatment, surgical resection with the maximum safe range in combination with postoperative concurrent radiotherapy and TMZ chemotherapy, the OS of patients with GBM remains unsatisfactory (25). Parsons et al. (26) found that IDH1 mutations occurred in a large number of young and secondary GBMs, and for the first time highlighted the potential of the mutant gene in the classification of GBM subtypes. A large number of in-depth studies on the biological behavior and molecular markers of GBM have improved our understanding of the pathogenesis of GBM. In 2016, the revised classification of central nervous system tumors was constructed based on genetic and epigenetic alterations, including IDH 1/2, MGMT methylation, 1p/19q codeletion, and EGFR (27). DNA methylation alterations, as a type of epigenetic mechanism, regulate gene expression levels via methylation status and are important in tumorigenesis and progression of GBMs (28,29). During the progression and recurrence of gliomas, DNA methylation is lost (30). The low methylation status leads to the increased expression of related oncogenes. GBM patients with MGMT promoter methylation are more sensitive to chemoradiotherapy and have excellent survival benefits (5). Klughamer et al. demonstrated that aberrant DNA methylation is associated with patient prognosis and can serve as a therapeutic target in GBM (31). In addition, considering the complexity and heterogeneity of GBMs, combined molecular markers are better than single biomarkers in the prognosis of glioma. Thus, DMDGs can be used as a tool for early diagnosis, risk stratification, and prognosis prediction, and can be therapeutic targets in GBM. Traditional DMDG screening occurs in both normal and tumor groups. In this study, screening events were conducted between low-grade astrocytomas (only astrocytomas were selected in LGG) and high-grade astrocytomas (GBM). Based on RNA-seq data of normal tissues in the GTEx database and gliomas in TCGA, 7006 significant DEGs were first selected between normal and glioma tissues. In the TCGA database, the 7006 significant DEGs were utilized to further screen for DEGs and AMGs in LGG and GBM. Seventy-two DMDGs were identified in GBM using the MethylMix algorithm. Function and pathway analyses showed that these genes were mainly enriched in the BP group and were mainly involved in cancerrelated pathways, including positive regulation of EMT, positive regulation of oligodendrocyte differentiation, and positive regulation of gliogenesis, suggesting that DMDGs play an (only astrocytoma) cohort (n = 195), and CGGA glioma cohort (n = 127). The AUC of the ROC curves indicated that the risk score had good prediction ability. Joint survival analysis combining the methylation values and the expression levels of DMDGs indicated that genes with high expression levels (MDK, SERPINB1, NMNAT33, PARVB, and UPP1) were significantly associated with poor prognosis and had the potential to be therapeutic targets in GBM. Among the six DMDGs, MDK (32,33), PDPN (34,35) and SERPINB1 (36) levels have been shown to be significantly higher in glioma tissues than in normal tissues. Enhanced levels of these genes is associated with the deterioration of prognosis in patients with GBM (37)(38)(39)(40). In addition, MDK (32,38,41,42), PDPN (43,44) and SERPINB1 (45) have been recognized as novel biomarkers and potent therapeutic targets for the treatment of GBM. In glioma, studies on NMNAT3, PARVB, and UPP1 are sparse. Only one bioinformatic analysis of UPP1 expression in gliomas has been published. It indicated that UPP1 expression is positively correlated with the grade of gliomas (46). Finally, because PARVB gene function in GBM has not yet been reported, it was further investigated. The mammalian parvin protein family includes three members (a-, b-, and g-Parvin) that have key roles in actin reorganization and focal adhesions (47). b-Parvin (PARVB) is overexpressed in human colorectal cancer (11) and tongue squamous cell carcinoma (12), and is closely associated with tumor progression. Importantly, overexpression of nuclear b-catenin and downregulation of Ecadherin were observed in human colorectal cancer featuring highly expression of PARVB protein (11). In PARVB knockdown oral cancer cells, loss of fibroblast-like structures at the wound edge was found (12). These observations are consistent with the potential close association between PARVB overexpression and EMT. Such an association between PARVB and EMT has not been previously reported in patients with GBM. EMT driven by signaling pathways is a biological process in which a polarized epithelial cell sheet undergoes multiple biochemical changes, ultimately leading to a mesenchymal phenotype characterized by cells with weak cell adhesion and enhanced migration (48). Moreover, cells undergoing EMT lose their apical basolateral polarization and acquire a fibroblast-like morphology, which increases their ability to migrate, spread, and disseminate to surrounding tumor tissues or distant sites (13,14,49,50). Function and pathway analyses showed that DMDGs (MDK, SERPINB1, NMNAT33, PARVB, and UPP1) were mainly involved in the positive regulation of EMT. Our study found that patients with higher PARVB tumor expression had significantly worse survival rates. The underlying mechanism by which PARVB promotes EMT was further investigated. GSEA of single genes demonstrated that highly expressed PARVB positively correlated with the JAK-STAT signaling pathway. Overexpression of PARVB in vitro can induce EMT, resulting in significantly increased cell proliferation, migration, and invasion. These findings imply that PARVB is involved in EMT-like processes in GBM. Silencing JAK2/STAT3 signaling partially inhibited the increased proliferation, migration, and invasion of PARVB-overexpressing cells in vitro. These data suggest that PARVB-induced EMT is at least partially mediated through the JAK2/STAT3 pathway.
Notably, this study differs in several ways from previous studies (20)(21)(22)(23). First, the study introduced a large number of normal groups (GTEx database) to overcome the limitation of analysis due to the lack of normal groups in TCGA. Second, it simplified the types of tumor subtypes. Only astrocytoma patients were chosen as the screening samples. This made the analysis more precise and targeted. Third, the search for DMDGs selected LGG and GBM, rather than between normal and tumor groups, which could amplify the small and easily overlooked differences in gliomas in previous studies. Fourth, the introduction of a correlation coefficient (Cor < -0.4) in our model more convincingly resulted in strong predictive power. Fifth, our DMDG prognostic model is suitable for GBM as well as being more suitable and accurate for LGG. Sixth, this is the first study to report PARVB in GBM. PARVB induced EMT, and affected cell proliferation, migration, and invasion through the JAK2/STAT3 pathway.

CONCLUSION
Combined analysis of transcriptome and DNA methylation profiles of TCGA was used to derive a 6-gene model-based risk score, screened from DMDGs. The approach allowed risk stratification survival prediction and personalized treatment plans for GBM. Our predictive model for DMDGs is better suited for LGGs than GBMs. Furthermore, we revealed that the PARVB gene can induce EMT to promote GBM cell malignant behavior partially through the JAK2/STAT3 pathway in GBM. Overall, we successfully constructed a prognostic model of DMDGS and validated the mechanism by which PARVB may be involved in regulating EMT in GBM.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Ethics Committee-approved study from the First Affiliated Hospital of Nanchang University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
All authors contributed to the experimental design and the analysis of data in this study. Conception and design: PW, BH. Acquisition, analysis and interpretation of data: WY, LM and FW. Writing, review, and/or revision of the manuscript: WY. Administrative, technical, or material support: BH, PW. Study supervision: ZJ, BH. All authors contributed to the article and approved the submitted version.