We collected 92 glioma tissues that were identified by the pathologists according to the 2016 WHO classification criteria from the Department of Pathology of the Affiliated Hospital of Xuzhou Medical University between 2016 and 2017. We obtained the ethical review and approval from the institutional ethics committee of Affiliated Hospital of Xuzhou Medical University (ethical review no. XYFY2018-KL056-01).

Human GBM cell lines (U118, U87, U251, T98G, and LN229), human normal brain glial cells (HEB), and normal human astrocyte (NHA) were originally obtained from the American Type Culture Collection (ATCC). All these cell lines were cultured with Dulbecco’s modified Eagle’s medium (DMEM; KeyGen Biotech, China) containing 10% fetal bovine serum (FBS, Takara, Japan) at 37°C within a humid atmosphere containing 5% CO₂.

We analyzed the MK3 expression with the “Single Gene Analysis” module of Gene Expression Profiling Interactive Analysis (GEPIA) web (http://gepia.cancer-pku.cn/) and observed the differential expression of MK3 between tumor and adjacent normal tissues for the different tumors or specific tumor subtypes of The Cancer Genome Atlas (TCGA). RNA-sequencing (RNA-seq) data of 689 GBM and low-grade glioma (LGG) tissues from TCGA dataset (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) and 1,157 normal tissues from Genotype-Tissue Expression (GTEx) (https://www.genome.gov/Funded-Programs-Projects/Genotype-Tissue-Expression-Project) were used for MK3 expression analyses. RNA-seq data of 413 glioma tissues (batch I) and 273 glioma tissues (batch II) datasets after deletion of incomplete data from the Chinese Glioma Genome Atlas (CGGA) dataset were also utilized for the expression analyses. In addition, we analyzed the expression of MK3 in glioma by using the Gene Expression Omnibus (GEO) database, GSE4290, and GSE7696 (https://www.ncbi.nlm.nih.gov/geo/). t-Test was used for the expression analyses by using R package “ggplot2.”

Cells were lysed with radioimmunoprecipitation assay (RIPA) extraction reagent [50 mM of Tris-HCl, pH 8.0, 150 mM of NaCl, 0.5% sodium deoxycholate, 1% NP-40, and 0.1% sodium dodecyl sulfate (SDS)] containing the protease inhibitors (Roche, Germany). Total protein was separated using 10% polyacrylamide gel electrophoresis and transferred to NC membrane (Millipore, USA). The membrane was blocked with 5% bovine serum albumin (BSA) in TBST for 1 h at room temperature and then incubated with the primary antibodies against MK3 (1:200, Santa Cruz, USA) and GAPDH (1:10,000, ProteinTech, USA) at 4°C overnight. Then the membrane was washed three times using 0.1% TBST buffer and incubated with the horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 h at room temperature, followed by three times TBST buffer washing. Finally, the signal was detected with enhanced chemiluminescence reagent.

The 4-μm-thick glioma sections were stained with a specific primary antibody against MK3 (1:50, Santa Cruz, USA) at 4°C overnight. Then the sections were incubated with secondary antibody (ZSGB-BIO, China) for 1 h and DAB for 2–5 min at room temperature, followed by hematoxylin staining. The images were acquired by Olympus microscopy and scored by two experienced pathologists without knowing the patients’ characteristics. Scores were calculated based on the intensity and percentage of positive tumor cells within the whole tissue, which were evaluated by using the German semiquantitative scoring method. The intensity score of cytoplasmic staining of each specimen was defined as follows: 0, negative; 1, weak; 2, moderate; and 3, strong. The percentage of positive cells was evaluated with 0.0% staining; 1, 1%–24 % staining; 2, 25%–49% staining; 3, 50%–74% staining; and 4, 75%–100% staining. The final immunoreactive score was calculated by the multiplication of the intensity scores and proportion scores.

We obtained the survival and clinical phenotype data of each sample from TCGA and CGGA datasets. The Kaplan–Meier (KM) curve analyses were conducted by the R package “survival” and “survminer.” Cox analysis was used for survival analyses by using the R package “survival.” The receiver operating characteristic (ROC) analysis was performed with the R package “pROC.”

Univariate and multivariate Cox analyses were conducted by using the R package “survival.” The WHO grade, IDH status, 1p/19q codeletion, primary therapy outcome, gender, age, and MK3 expression were included in these analyses.

We searched the MK3 with the query of protein name (“MK3”) and organism (“Homo sapiens”) on the STRING website (https://string-db.org/). The basic setting parameters were the network type (“full STRING network”), meaning of network edges (“evidence”), active interaction sources (“experiments”), minimum required interaction score (“low confidence (0.150)”), and max number of interactors to show (“no more than 50 interactors in 1^st shell”).

Gene expression data of GBM and LGG in HTSeq-Counts were downloaded from TCGA website for further analysis. The correlated genes of MK3 were screened with Pearson’s correlation coefficients (|r| >0.4 and p < 0.001) with R package “deseq2.” Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genome (KEGG) analysis were conducted on MK3 correlated genes with R package “clusterProfiler” to identify the possible biological functions and signaling pathways affected by MK3. The biological process (BP), cellular component (CC), and molecular function (MF) were applied in GO analysis.

The Gene Set Enrichment Analysis (GSEA) was performed to dissect the cancer-related pathways by using the MSigDB Collection (c2.cp.v7.2.symbols.gmt) of the clusterProfiler R package. The gene sets with normal p. adjust value and false discovery rate (FDR) q value both less than 0.05 were considered as significantly enriched.

We obtained the single-cell RNA-seq data from the Single Cell^BETA PORTAL website (https://singlecell.broadinstitute.org/single_cell/study/SCP393/single-cell-rna-seq-of-adult-and-pediatric-glioblastoma#study-download), (https://singlecell.broadinstitute.org/single_cell/study/SCP50/single-cell-rna-seq-analysis-of-astrocytoma), and (https://singlecell.broadinstitute.org/single_cell/study/SCP147/single-cell-analysis-in-pediatric-midline-gliomas-with-histone-h3k27m-mutation#study-visualize) (26) to explore the distribution of MK3 from the tumor samples.

The immune cell infiltration levels from tumors were quantified by the single-sample GSEA (ssGSEA) with the R package “GSVA” (27). The correlation between the MK3 expression and immunoinhibitors, chemokines, and chemokine receptors across human cancers was determined on the Tumor-Immune System Interaction Database (TISIDB) website (http://cis.hku.hk/TISIDB/index.php). The correlation analysis between the expression of MK3 and immunoinhibitors, chemokines, and chemokine receptor genes of the glioma was evaluated by using the R-package “ggplot2.”

The correlation of MK3 expression with immune checkpoint expression from glioma was analyzed on the Tumor IMmune Estimation Resource (TIMER) website (https://cistrome.shinyapps.io/timer/) and adjusted by the tumor purity.

All gene expression data were normalized with log2 transformation. Statistical analyses were conducted by R (3.6.3). The survival data from the CGGA database were acquired by KM. The correlation analysis was evaluated by chi-square (χ²) test, Pearson’s correlation, or Spearman’s correlation analysis. Student’s t-test was used to determine the statistical significance of differences between groups, and the differences between more than two groups were analyzed by one-way ANOVA or the Kruskal–Wallis test. A two-sided Mann–Whitney test was used for non-parametric data. p-Value < 0.05 was considered statistically significant. *p < 0.05, **p < 0.01, and ***p < 0.001.

To determine the MK3 expression in glioma, we analyzed the GEPIA database that contains the gene expression profile information across all tumor samples and paired normal tissues. The MK3 was aberrantly overexpressed in GBM (n = 163) and LGG (n = 518) tissues compared with the normal brain tissues (n = 207) that had a low level of MK3 mRNA expression among different normal human organs ( Figures 1A, B ). We further validated the data by using the GEO database; a similar overexpression pattern was observed in the glioma tissues compared with the non-tumor (GSE4290) or normal (GSE7696) tissues ( Figures 1C, D ). We also assessed the MK3 expression levels in glioma tissues by immunohistochemistry (IHC) staining. The staining results showed a similar MK3 expression pattern from the tumor tissues when compared with the para-tumor tissues ( Figures 1E, F ).

Furthermore, we measured the protein expression levels of MK3 in five glioma cell lines (U118, U87, U251, T98G, and LN229), HEB, and NHA cell lines by Western blotting analysis. We observed that the protein levels of MK3 were also higher in the majority of the glioma cell lines ( Figure 1H ). Collectively, our results suggest that MK3 is aberrantly overexpressed in glioma.

The significant positive correlation between the MK3 level and WHO grade in glioma tissues ( Figures 1E, G , Table 1 ) suggested a potential role of MK3 in the clinicopathologic features of glioma patients. To determine the clinical implication of MK3 expression, we segmented the patients of GBM and LGG in TCGA database into low or high MK3 expression groups according to the median value. As shown in Table 2 , MK3 overexpression was significantly correlated with the WHO grade, IDH status, 1p/19q codeletion, and age. We also analyzed the expression of MK3 with OS, progression-free interval (PFI), and disease-specific survival (DSS) event of GBM and LGG, and the results showed that the expression of MK3 was lower from alive patients ( Figures 2A–C ). Moreover, we determined the clinical implication of MK3 expression from human GBM and LGG samples in TCGA database, and the significant differences were found in 1p/19q codeletion ( Figure 2D ), WHO grade ( Figure 2E ), and IDH status ( Figure 2F ). The results were also validated with the CGGA database ( Figures 2G–L ). Altogether, we demonstrate that MK3 expression is closely correlated with the poor clinicopathologic features of glioma patients.

To investigate the relationship between the MK3 expression and the prognosis of glioma patients, we divided the cancer cases into high-risk and low-risk groups according to the cutoff value of the median risk score. The poorer prognosis and higher death rate were observed in the high-risk group ( Figure 3A ). Glioma patients were also separated into the MK3-high and MK3-low groups based on the expression levels with the median value to generate a KM survival curve. The KM survival curve demonstrated that the OS rate of glioma patients in the high-expression group was significantly poorer than that in the low-expression group in both the CGGA and TCGA datasets ( Figures 3B–D ). The PFI and DSS rate of glioma patients were also significantly and negatively correlated with MK3 expression in TCGA dataset ( Figures 3E, F ). In addition, we performed a subgroup survival analysis of OS with TCGA dataset, and we found that high expression of MK3 was associated with poor prognosis in WHO grade G3, non-codel 1p/19q codeletion, age less than or equal to 60 years, age greater than 60 years, and histological type astrocytoma subgroup of glioma ( Figures 3G–O ). We also conducted the ROC curve analysis to evaluate the diagnostic value of MK3, and the area under the curve (AUC) was 0.966, which indicated a high diagnostic value of MK3 in glioma ( Figure 3R ).

To further evaluate the prognostic value of MK3, we performed the univariate and multivariate Cox regression analyses, which are the widely used approaches for identifying predictive biomarkers. Univariate Cox analysis unearthed that MK3 was a high-risk factor (hazard ratio (95% CI) = 3.037 (2.349–3.297), p < 0.001). Moreover, WHO grade, IDH status, 1p/19q codeletion, primary therapy outcome, age, and MK3 expression were all significantly correlated with poor OS ( Table 3 ), which indicated their association with the OS in glioma. On the other hand, the multivariate Cox analysis demonstrated that only WHO grade, IDH status, primary therapy outcome, gender, and age were independent prognostic factors for glioma ( Table 3 ). In conclusion, these results demonstrate that MK3 might serve as a valuable prognostic biomarker for glioma patients.

To identify the functional and physical interaction partners of MK3, we used STRING to analyze the protein–protein interactions. As shown in the interaction network ( Figure 4A ), a total of 38 MK3 interaction proteins were found by experimental evidence. Furthermore, we recovered 3,067 positively correlated genes (Pearson’s correlation coefficient r > 0.4) and 1,229 negatively correlated genes (Pearson’s correlation coefficient r < −0.4) of MK3 from TCGA transcriptome dataset. The top 15 positively correlated genes and negatively correlated genes of MK3 were plotted in the heatmap ( Figure 4B ). To further uncover the potential role of MK3 in glioma, we performed GO and KEGG analyses by TCGA dataset on MK3 correlated genes. Interestingly, the results of GO classification demonstrated that MK3 correlated genes were significantly enriched in immune-related functions among the top-ranked BPs, including the neutrophil activation, neutrophil-mediated immunity, neutrophil activation involved immune response, T-cell activation, regulation of innate immune response, leukocyte cell–cell adhesion, regulation of T-cell activation, leukocyte proliferation, antigen processing and presentation, response to interferon gamma, antigen processing and presentation of peptide antigen, cellular response to interferon gamma, T-cell proliferation, and interferon gamma-mediated signaling pathway ( Figures 4C, D ). KEGG pathway analysis showed that MK3 correlated genes were also closely associated with immune-related signaling pathways, including cytokine–cytokine receptor interaction, osteoclast differentiation, chemokine signaling pathway, proteoglycans in cancer, Th17 cell differentiation, natural killer cell-mediated cytotoxicity, TNF signaling pathway, Toll-like receptor signaling pathway, Th1 and Th2 cell differentiation, T-cell receptor signaling pathway, B-cell receptor signaling pathway, and NF-kappa B signaling pathway ( Figures 4E, F ).

We also performed GSEA by using the MK3-low and MK3-high datasets from TCGA transcriptome dataset to identify signaling pathways that are affected by MK3 overexpression in glioma. Gene sets related to chemokine signaling, Toll-like receptor signaling, cytokine–cytokine receptor interaction, natural killer cell-mediated cytotoxicity, interferon signaling, interferon gamma-signaling, cell surface interactions at the vascular wall, neutrophil degranulation, extracellular matrix organization, signaling by interleukins, signaling by the B-cell receptor, and Toll-like receptor cascades pathways were significantly enriched in the high MK3 expression group ( Figures 5A–L ). Together, our results suggest that MK3 might participate in the tumor immune microenvironment and immune regulation.

Considering the tumor heterogeneity, we further analyzed the distribution of MK3 expression in glioma samples by assessing single-cell RNA-seq data. We observed that the extensive expression of MK3 was in immune cells (macrophage and T cell) ( Figures 6A–C ), which further support our IHC staining data in glioma tissues ( Figure 1E ). These results suggest the potential role of MK3 in regulating immune response. Therefore, we determined the correlation between MK3 expression and immune infiltration in glioma by ssGSEA with Spearman’s R. We found that MK3 expression was negatively correlated with infiltrating levels of T follicular helper (TFH), CD8 T cells, T gamma delta (Tgd), T central memory (Tcm), NK CD56bright cells, and plasmacytoid DC (pDC) cells and was positively correlated with infiltrating levels of macrophages, neutrophils, eosinophils, activated dendritic cells (aDC), immature DC (iDC), T, Th17, cytotoxic cells, NK CD56dim, Th2, NK, T helper, Th1, B cells, and DC (all p < 0.05, Figures 6D, E ). These data indicate that MK3 might participate in immune infiltration in glioma.

Taking advantage of TISIDB website, we also explored the correlation between the expression of MK3 and immune-related genes in glioma. Our results suggest that the expression of MK3 was positively correlated to the majority of immunoinhibitors ( Figures 7A, B ), pro-tumor chemokines ( Figures 7C, D ), and chemokine receptors ( Figures 7E, F ) in both GBM and LGG patients.

Programmed cell death protein 1 (PD-1/PDCD1), programmed cell death 1 ligand 1 (PD-L1/CD274), and cytotoxic T lymphocyte antigen 4 (CTLA-4) are vital immune checkpoints that play important roles in tumor immune escape. And they also served as predictive markers for the therapeutic efficacy of immune checkpoint inhibitors (ICIs). Thus, we especially assessed the relationship of MK3 with PD-1, PD-L1, and CTLA-4 in glioma by using TIMER tool. The results from TIMER website illustrated that there was a significant positive correlation between MK3 expression and PD-1, PD-L1, and CTLA-4 in glioma, which was adjusted by tumor purity ( Figures 8A–F ). Similarly, we also found that the expression of MK3 was significantly positively correlated with PD-1 (r = 0.520, p < 0.001), PD-L1 (r = 0.570, p < 0.001), and CTLA-4 (r = 0.360, p < 0.001) in glioma ( Figures 8G–I ). Taken together, these results suggest that MK3 might participate in tumorigenesis and the development of glioma by regulating the tumor immune escape.