A total of 184 surgically resected human glioma specimens were collected in the Department of Neurosurgery, the First Affiliated Hospital, Zhengzhou University, between 2011 and 2019. Among 184 patients, 115 (62.5%) were males and 69 (37.5%) were females. According to the WHO classification of CNS tumors (Wesseling and Capper, 2018), the collected neoplasms were classified as follows: 144 patients with grade III gliomas (including 2 astrocytoma, 109 anaplastic astrocytoma, 17 anaplastic oligoglioma, 4 oligoastrocytoma, 3 anaplastic oligoastrocytoma, and 9 anaplastic ependymoma) and 40 patients with GBM (grade IV). Collected clinical data included the gender, tumor location, histology, cell origin, tumor grade, clinical surgery status, and radiochemotherapy status. All tissues were snap-frozen in liquid nitrogen and stored at −80°C after resection. The clinicopathological features and surgery status of the patients are summarized in Table 1.

To block the endogenous peroxidase activity, the sections from 4% formalin-fixed and paraffin-embedded tissue blocks were dewaxed in xylene, rehydrated using graded concentration ethanol, and soaked in 0.3% hydrogen peroxide, which were cut into 5 μm-thick sections. Then, the tissue sections were rinsed in PBS. After having blocked in 10% goat serum (Gibco, Waltham, MA, USA) for nonspecific reactions, the sections were incubated with the anti-PLAGL2 polyclonal antibody (1:100, Abcam, Cambridge, MA, USA) overnight at 4°C. Negative control slides were processed in parallel with a nonspecific IgG. After washing, the sections were incubated with horse reddish peroxidase (HRP)-labeled anti-rabbit secondary antibody (Santa Cruz, CA, USA) for 1 h at room temperature. The slides were subjected to staining with 3, 30-diaminobenzidine solution and briefly counterstained with hematoxylin. Finally, the sections were rehydrated with xylene and ethanol, mounted in mounting resin, and examined with a microscope.

Two independent pathologists blinded to clinicopathological outcome calculated immunostaining reactions by multiplying the intensity and proportion of positive tumor cells. Briefly, the scores about the proportion of positive cells (0, <5%; 1, 6–50%; 2, >50%) and the scores about the staining intensity (0, no or weak staining; 1, moderate staining; and 2, strong staining) were combined. All the cases were divided into low expression and high expression subjected to the total score (median). Any discrepancy in the scoring was resolved by discussion between two pathologists in less than 10% of the examined slides, and an agreement was reached by rescoring the sections combined with a third pathologist in a multi-viewer microscope.

Cancer Cell Line Encyclopedia (CCLE; https://portals.broadinstitute.org/ccle) was downloaded from Gene Expression Omnibus (GEO, series GSE36133) and was used to identify the alternation of the expression of PLAGL2 across various cancer types. The Affymetrix Human Genome U133 Plus 2.0 DNA microarray gene expressions of 84 CNS cell lines were downloaded from CCLE in March 2012. Robust multiarray average (RMA) normalization and gene expression values was performed. One-way ANOVA was used to compare the expression of PLAGL2 in different types of cancer cell lines (Barretina et al., 2012; Cancer Cell Line Encyclopedia Consortium and Genomics of Drug Sensitivity in Cancer Consortium, 2015).

The transcript level of PLAGL2 in glioma was ascertained by the ONCOMINE database (https://www.oncomine.org/), with a threshold set as such that P < 1E − 4, fold change >2, top gene rank 10% (Rhodes et al., 2004). The mRNA levels of PLAGL2 in cancer tissues were in contrast with that in normal tissues. Two one-sided T-tests were used to evaluate the differences. The cut-off of p-value was defined as 0.01, and the cut-off of fold change value was identified as 1. The microarray data of patients with glioma were downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo) public database under accession numbers GSE4290 and GSE53733 (Sun et al., 2006; Reifenberger et al., 2014).

Normalized RNAseq expression data and corresponding clinical material were obtained from two independent databases: The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) were both downloaded from Gliovis (http://gliovis.bioinfo.cnio.es/) (Bowman et al., 2017). TCGA is a public database (http://cancergenome.nih.gov/), which includes 29 cancer types, along with gene expression data and clinical information. Furthermore, a database (ID: mRNAseq 325) consisting of 325 glioma tissues with different grades (WHO I–IV) was downloaded from the CGGA (Wang et al., 2015; Liu et al., 2018). LGG and HGG were defined as WHO grades I–II and WHO grades III–IV according to the 2016 WHO classification of CNS tumors, respectively (Louis et al., 2021).

The UALCAN (http://ualcan.path.uab.edu) and GEPIA datasets (http://gepia.cancer-pku.cn/) of TCGA gene expression were used to analyze the expression of PLAGL2 in GBM tissues and normal tissues (Asati et al., 2016). Two one-sided t-tests were used to evaluate the differences.

PLAGL2 were analyzed at the protein and RNA levels in the Human Protein Atlas (HPA) dataset (https://www.proteinatlas.org/), which contained the immunohistochemistry staining images of 153 clinical GBM patients and their survival information with 54 females and 99 males (Pontén et al., 2008; Uhlen et al., 2017).

The Tumor Immune Estimation Resource (TIMER2.0) (http://timer.cistrome.org/) is a public web database for comprehensive assessments of abundance of tumor-infiltrating immune cells (TIICs) according to gene expression profiles (Li et al., 2020). We mainly performed PLAGL2 expression in GBM and the correlation between gene expression and abundance of immune infiltrates, involving B cells, CD8⁺ T cells, and CD4⁺ T cells, macrophages, neutrophils, dendritic cells (DCs), compared with purity in GBM and assessed how PLAGL2 is correlated with immune cell markers including CD8+T cells, T cells, B cells, monocytes, M1 macrophages, M2 macrophages, neutrophils, DCs, Th1 cells, type 2 helper T cells (Th2), Tfh cells, type 17 helper T cells (Th17), and Treg.

GENT2 analysis (http://gent2.appex.kr/gent2/) revealed the gene expression in several types of cancer. The data from about 49,000 healthy and cancer individuals were analyzed in more than 30 different cancer types (Park et al., 2019).

The GeneMANIA (http://www.genemania.org/) web database revealed the interaction relationship of a gene and the performed gene function, physical interaction, pathway, co-expression, and localization (Warde-Farley et al., 2010). The STRING database (https://string-db.org/) was utilized to obtain the interaction network of related genes for PLAGL2, including direct (physical) and indirect (functional) associations. Thus, nodes represent genes, whereas links conduct connected networks (Szklarczyk et al., 2017; Szklarczyk et al., 2019).

To investigate the potential mechanisms underlying the interaction of PLAGL2 expression on glioma progression, a GSEA was conducted to screen out whether some biological pathways showed statistically significant differences between high and low PLAGL2 expression groups (Subramanian et al., 2005). Pre-defined gene sets were obtained from the Molecular Signatures database, MSigDB (http://software.broadinstitute.org/gsea/msigdb) (Subramanian et al., 2005). Gene sets: M457 (CYTOSKELETON), M2001 (WU_CELL_ MIGRATION), and M16210 (CELL_PROLIFERATION_GO_0008283), M14493 (Angiogenesis_GO_ 0001525), M39729 (WP_VEGFAVEGFR2_ SIGNALING_PATHWAY) (Liberzon et al., 2011). Normalized enrichment score (NES) and false discovery rate (FDR) were calculated to verify the significant difference for GSEA. For each analysis, gene set permutations were implemented for 1,000 times. Gene sets with an FDR <0.05 and normal p <0.05 were viewed as significantly enriched.

Statistical analyses were performed by SPSS 22.0 software, and the results were presented as the mean ± standard deviation (SD). Student’s t-test was performed to detect significant differences between groups. The relationships of PLAGL2 expression and various clinicopathological features were estimated by chi-squared test. Analyses with more than two groups were performed by the one-way ANOVA test. The Kaplan–Meier with log-rank test was used to analyze the survival data. The correlation analysis between PLAGL2 expression and immune-related markers was performed using Pearson correlation. Univariate and multivariate Cox regression analyses were used to analyze the clinical prognostic parameters and related independent risk factors. The results were considered significant at *p < 0.05, **p < 0.01, and ***p < 0.001.

We performed pan-cancer analysis derived from TCGA by adopting Timer 2.0 and Gent2 so as to estimate the expression pattern of PLAGL2. According to Timer analysis, PLAGL2 was upregulated in bladder, breast, cervix, bile ducts, colorectal, esophagus, kidney, liver, lung, stomach, and CNS cancer compared with normal cases (Figure 1A, Supplementary Table S1), whereas it was downregulated in renal and thyroid cancers.

Moreover, we estimated the PLAGL2 expression data profiles of 72 paired cancers vs. normal tissues by utilizing the HG-U133 microarray (GPL570 platform) of the Gent2 database. The results showed that PLAGL2 was upregulated in several cancers, including adrenal gland, brain, breast, cervix, colon, esophagus, liver, lung, oral, ovary, and skin. Conversely, the upregulation of PLAGL2 was reported in the pharynx, small intestine, testes, stomach, head, and neck (Figure 1B, Supplementary Table S2). Combining the above two results, systematic mRNA expression analysis with the Timer and GENT2 databases have shown that PLAGL2 was consistently upregulated across a wide range of cancer types that included breast, cervix, colon, esophagus, liver, lung, and brain.

Next, we used these respective cancer types in which significant upregulation was found for further investigation. In addition, we selected CCLE and HPA databases to further validate the PLAGL2 expression in the cell line. The CCLE analysis results indicated a certain degree of PLAGL2 expression in 84 CNS cell lines containing glioma cancer cell lines (Supplementary Table S3). In GBM U251-MG cells from HPA, the immunofluorescent staining of human cell line U251-MG revealed that PLAGL2 was localized to the nucleoplasm and cytosol, as shown in Supplementary Figure S1, where green represents antibody, red means microtubules, and blue means nucleus.

Firstly, we assessed any deregulated differences of the PLAGL2 mRNA level in glioma compared to normal brain tissues, which were tested using the data in human glioma from publicly available microarray datasets (https://www.oncomine.org/) (Rhodes et al., 2004). Based on the ONCOMINE database, PLAGL2 has significantly high expression in glioma tissues (astrocytoma and glioblastoma) as compared with normal tissues (Figures 2A,B) (Rickman et al., 2001; Murat et al., 2008). Similarly, as shown in Figure 2C and Supplementary Figures S2A,B, PLAGL2 are upregulated in GBM from the TCGA sample based on UALCAN and GEPIA databases. Meanwhile, the results in LGG indicated that patients with glioma also had higher expression levels of PLAGL2 compared with non-tumor tissue, which indicated that PLAGL2 expression was closely linked to the malignancy of glioma (Supplementary Figure S2C). To further characterize the relationship of PLAGL2 expression level and tumor specimens, we explored the expression levels of PLAGL2 in the GEO database (GSE4290). According to the WHO grade, the results showed that the expression of PLAGL2 mRNA was greatly higher in grades III and IV than in normal brain tissues (Figure 2D). Moreover, GENT2 database analysis showed that PLAGL2 expression increases with tumor grade (Supplementary Figure S3).

IDH mutation and 1p/19q co-deletion status are well-known clinically relevant molecular markers in glioma. Herein, we analyzed the correlation between the mRNA levels of PLAGL2 and major clinical features of GBM patients in CGGA, including grade level, IDH mutation status, 1p/19q co-deletion status, gender, and age. As shown in Supplementary Figures S4A−E, PLAGL2 were differently expressed in different grade groups and significantly associated with IDH mutation, 1p/19q co-deletion, gender, and age status, respectively. Overall, the above results illustrated that the mRNA level of PLAGL2 was significantly associated with clinical and pathological features, such as grade, gender, age, IDH mutation, and 1p/19q co-deletion status, which were also considered as potential prognostic and diagnostic predictors of HGG with poorer treatment.

As described above, the mRNA levels of PLAGL2 were upregulated in GBM, so we hypothesized that the protein levels of the PLAGL2 were also elevated. In the HPA dataset, the IHC staining data revealed that PLAGL2 was mainly localized in the nucleus and had higher expressions in GBM tissues (Figure 3A). To validate the above findings and investigate the clinicopathological roles and distribution of PLAGL2 expression in HGG, an immunohistochemical analysis of the 184 paraffin-embedded HGG tissue blocks was performed. These gliomas comprised 144 grade III and 40 grade IV tumors. Representative immune-histochemical staining of PLAGL2 in gliomas is illustrated in Figure 3B (a-d). PLAGL2 expression levels were significantly lower in grade III (a-b) compared to grade IV gliomas (c-d). A strong expression of PLAGL2 was observed in the nuclei of tumor cells. Additionally, there was no positive signal in any non-neoplastic tissues. The cases were divided into high or low PLAGL2 expression group by the median of PLAGL2 immunostaining score in cancer tissues (n = 92 for each group). As listed in Table 1, the expression of PLAGL2 was significantly correlated with WHO grade (p < 0.05), but not with gender, tumor locus, histology, cell origin, treatment regimen, and radiochemotherapy. Meanwhile, considering the GBM group, the analysis results showed that the expression of PLAGL2 was significantly correlated with gender (p < 0.01), but not with tumor locus and radiochemotherapy (Table 2). Those results suggested that PLAGL2 expression was more prevalent in aggressive gliomas.

Survival length was determined from the time of primary tumor surgery to the time of death or the last follow-up. To investigate the effect of PLAGL2 overexpression on the overall survival (OS) and progression-free survival (PFS) of HGG patients and explore the prognostic value of PLAGL2 in HGG, survival analysis was performed based on information of 184 HGG patients. As shown in Figures 4A,B, HGG patients with high PLAGL2 expression showed significantly worse OS and PFS than those with low PLAGL2 expression (p < 0.001 and p = 0.006, respectively). Similarly, we explored the expression levels of PLAGL2 in the GEO database. GSE53733 experiments indicated that the expression of PLAGL2 in the short-term OS group was significantly higher than that in the long-term OS group (Figure 4C). Moreover, we used the Gent2 database to analyze the survival plot or KM plot for PLAGL2 expression against HGG (Figure 4D). The GENT2 database was used to obtain the log-rank test curve, which showed that high levels of the PLAGL2 expression group had a lower survival rate in HGG patients (p < 0.05). Moreover, a low prognostic value was observed in high-grade stage groups (p = 0.001) (Figure 4E). Next, HGG patients with primary glioma and recurrent glioma in CGGA were grouped into cohorts based on the median expression levels of PLAGL2 expression. Similarly, we found that those with low expression had a better prognosis compared to those with higher expression levels in primary and recurrent glioma (Supplementary Figures S5A, 5B). A higher expression level of PLAGL2 also predicted a poor prognosis for WHO grade III and IV glioma patients compared to those with lower expression (Supplementary Figures S5C−F). Furthermore, the AUC values of the prognostic model for the 1-, 2-, and 3-year survival rate prediction in the TCGA cohort were 0.708, 0.715, and 0.713, respectively (Supplementary Figure S6A). The AUC values for disease specific survival in the TCGA cohort at 1, 2, and 3 years were 0.713, 0.72, and 0.714, respectively (Supplementary Figure S6B). These findings confirmed that the prediction performance of the prognostic signature was considerably promising and the survival analysis accurately predicted the prognosis of patients with GBM and LGG. An ROC curve analysis was also performed in that the high sensitivity and specificity of differentiating IDH status, histological, and 1p19q co-deletion status are based on the OS of the TCGA database (Supplementary Figure S7).

To further determine the prognostic value of PLAGL2 expression, we performed the univariate and multivariate Cox regression analyses based on OS and PFS of HGG patients. As shown in Tables 3, 4, PLAGL2 expression, age class, and WHO grade resulted as independent prognostic factors for the OS and PFS of HGG patients. Next, the results of univariate and multivariate Cox regression analyses suggested that clinical histology, age, and 1p19q co-deletion status were associated with the OS of HGG patients ( Supplementary Table S4 ).

Tumor-infiltrating lymphocytes (TILs) could affect a patient’s OS and regulate the tumor response to therapies (Park et al., 2020). By combining the above results, we found that PLAGL2 was significantly upregulated and associated with OS and PFS in HGG patients. Cumulative studies have suggested that PLAGL2 expression has a vital role in the malignancy of various cancers (Zheng et al., 2010; Zhao et al., 2020; Wu et al., 2021). To explore the pan-cancer correlation between PLAGL2 expression and immune infiltration, we first evaluated the abundance of immune cell infiltration. As shown in Figure 5, we adopted the TIMER database to illustrate the profiles of PLAGL2 correlating with various immune infiltrations, which showed that it was positively correlated with the immune cell infiltration levels of neutrophils. However, data indicated that it was obviously contradicted in B cells, macrophages, CD4⁺T cells, CD8⁺T cells, and DCs. Next, we then performed a correlation analysis between PLAGL2 and tumor-infiltrating immune cells, including B cells, CD4⁺T cells, CD8⁺T cells, macrophages, neutrophils, and DCs by the TIMER database in GBM patients to evaluate the immunotherapy effect. Tumor purity acted as a critical factor, which affected the analysis of immune infiltration in the analysis of the genomic approach. As shown in Figure 6A, PLAGL2 was positively associated with purity (rho = 0.265, p = 0.0017), B cells (rho = 0.178, p = 0.0371), CD4⁺T cells (rho = 0.335, p < 0.001), macrophages (rho = 0.231, p < 0.01), neutrophils (rho = 0.461, p < 0.001) and DCs (rho = 0.428, p < 0.001), and negatively associated with CD8⁺T cells (rho = -0.314, p < 0.001) in GBM. Additionally, we also investigated the immune infiltration status of PLAGL2 related to immune cells in LGG. Similarly, the purity (rho = 0.265, p = 0.0125), B cells (rho = 0.148, p = 0.00115), CD4⁺T cells (rho = 0.253, p < 0.001), CD8⁺T cells (rho = -0.14, p = 0.00218), macrophages (rho = 0.136, p = 0.00281), neutrophils (rho = 0.495, p < 0.001), and DCs (rho = 0.271, p < 0.001) were significantly correlated with PLAGL2 in LGG (Figure 6B). In general, our findings revealed distant associations between PLAGL2 expression with cancer purity and immune cells infiltration, such as B cells, CD8⁺T cells, CD4⁺T cells, and macrophages and neutrophils, while DCs revealed that PLAGL2 had a vital role in immune infiltration in GBM, which modulated the infiltration of immune cells to tumor tissues.

Next, we comprehensively clarified the correlation of PLAGL2 and the status of infiltrating immune cells in GBM based on a cohort of immunological makers using the TIMER2.0 database. The gene markers of immune cells were used to analyze and identify the immune cells, including CD8⁺ T cells, CD4⁺ T cells, regulatory T cells, B cells, tumor-associated macrophages (TAMs), M1 and M2 macrophages, monocytes, neutrophils, DCs, and natural killer (NK) cells. Meanwhile, NK cells, T helper 1 (Th1), T helper 2 (Th2), follicular helper T (Tfh), T helper 17 (Th17), and regulatory T (Tregs) were also analyzed. The correlation analysis data were adjusted for tumor purity. An analysis of the TIMER database showed that PLAGL2 expression in GBM was significantly correlated with the expression of marker genes in CD8⁺T cells, general T cells, B cells, monocytes, M1 and M2 macrophages, neutrophils, DCs, NK cells, Th1, Th2, Tfh, Th17, and Treg (Table 5). Interestingly, PLAGL2 expression was critically correlated with the expression of markers of specific immune cells such as CD8⁺T cell marker, CD8B (r = -0.314, p < 0.001); T-cell marker, CD3D (r = -0.284, p < 0.001) and CD2 (r = -0.185, p < 0.05); B-cell marker, CD79B (r = -0.178, p < 0.05); monocyte marker, CD14 (r = 0.171, p < 0.05) and CD115 (r = 0.175, p < 0.05); M1 macrophage marker, CCL5 (r = -0.184, p < 0.05), IL18 (r = -0.266, p < 0.01) and COX (r = 0.396, p < 0.001); M2 macrophage marker, CD163 (r = 0.178, p < 0.05), neutrophil marker, CD11b (r = 0.337, p < 0.001); DC marker, HLA-DQB1 (r = 0.221, p < 0.01) and BDCA4 (r = 0.374, p < 0.001), and NK cell marker, NCR3 (r = -0.219, p < 0.01).

Moreover, the expression of PLAGL2 was associated with the expression of markers of specific subsets of T cells in GBM, which included the Th1 marker, TBX2 (r = 0.407, p < 0.001); Th2 marker, CD14 (r = 0.171, p < 0.05), GATA3 (r = 0.177, p < 0.05), STAT6 (r = 0.265, p < 0.01), STAT5A (r = 0.43; p < 0.001); Tfh marker, BCL6 (r = 0.444, p < 0.001); Th17 marker, STAT3 (r = 0.608, p < 0.001); Treg marker, FOXP3 (r = 0.231, p < 0.01), STAT5B (r = 0.612, p < 0.001), and TGFB (r = 0.363, p < 0.001). However, PLAGL2 expression did not reveal any significant correlation with CD4⁺ T cells and TAMs. In summary, a comprehensive analysis indicated that PLAGL2 expression was strongly correlated with infiltration of immune cells in GBM.

An aberrant expression of PLAGL2 may be involved in the process of diverse cancer types. Herein, we implemented two web-based network tools, GeneMANIA, and STRING, to investigate the interaction network related to PLAGL2. It is well known that protein mediates a wide variety of cellular functions and biological processes regarding protein interaction and the signal transduction pathway (Pawson and Nash, 2000; Westermarck et al., 2013). GeneMANIA, as an integrated network, focuses on functional prediction and performs an interaction network analysis based on a network-based gene ranking algorithm. Meanwhile, STRING is more inclined to the physical and functional interactions of a gene set. As Figure 7A showed, GeneMANIA supplied that the PPI-associated protein included ORC6, CMTR1, KIF3B, LMNB2, PGD, NEK2, MTHFD2, PTTG1, EZH2, NCAPH, PKMYT1, RCHY1, SRPK1, RBM38, DMRTC2, ASXL1, SF1, PLAS2, PIAS4, and PIAS1. STRING analysis provided a predicted PPI network, which shares interaction with RBBP8NL, PDPN, TM9SF4, SOX15, POFUT1, KIF3B, RCHY1, ENC1, RAG1, and GPR85 (Figure 7B). The PPI network stats are the number of nodes: 11; the number of edges: 17; average node degree: 3.09; local clustering coefficient: 0.893, and PPI enrichment p < 0.05. The related parameters predicted that PLAGL2 is involved in the progression and prognosis of cancer.

To investigate the biological characteristics shared by the different PLAGL2 expression levels on proliferation, migration, and F-actin polymerization of gliomas, we divided HGG patients from GSE4290 into PLAGL2-positive and PLAGL2-negative groups and performed GSEA, a robust computational method that determines whether an a priori defined set of genes is statistically significant, as well as the concordant differences between both groups. The enrichment plots of GSEA showed that the gene signatures of proliferation were enriched in PLAGL2-positive glioma HGG tissues compared with PLAGL2-negative groups (Figure 8A). In addition, similar enrichment was found related to the gene signatures of migration when the PLAGL2-positive glioma HGG tissues were compared with PLAGL2-negative groups (Figure 8B). Furthermore, the cytoskeleton assembly gene sets were enriched in PLAGL2-positive glioma HGG tissues as compared with PLAGL2-negative groups (Figure 8C). Accordingly, GSEA results further indicated that glioma with various PLAGL2 mRNA expression levels had the distinct status of angiogenesis (Figure 8D). The VEGF-VEGFR2 signaling pathway represented a growth factor with important pro-angiogenic activity, having a mitogenic and anti-apoptotic effect on endothelial cells, increasing the vascular permeability, and promoting angiogenesis (Ferrara et al., 2003). The gene signatures of the VEGFA-VEGFR2 signaling pathway assembly were abundantly activated in patients with a higher expression of PLAGL2 (Figure 8E).