The RNA-seq data and corresponding clinical information were downloaded from CGGA (www.cgga.org.cn) (38), TCGA (http://xena.ucsc.edu/) and Rembrandt (https://wiki.cancerimagingarchive.net/display/Public/REMBRANDT) database, including 713 samples from CGGA database (20 normal samples and 693 glioma samples), 1680 samples from TCGA database (895 normal samples and 685 glioma samples) and 364 samples from Rembrandt database (21 normal samples and 343 glioma samples). Clinical information of the glioma patients consisted of WHO grade, IDH1 status, 1p/19q status, age, gender and overall survival. Some samples with unavailable or unclear clinical information were removed. In addition, the DNA methylation data was also downloaded from TCGA database described above.

All glioma samples were divided into high and low MD2 expression (high methylation or low methylation) groups by the median expression level of MD2 in each database. The association between MD2 expression level (or methylation level) and overall survival in glioma samples were assessed by Kaplan-Meier survival analysis with log-rank test.

Univariate and multivariate Cox regression was used to examine whether MD2 expression, age, pathological grade, 1p/19q status and IDH mutation were independent prognostic factors in glioma patients based on CGGA and TCGA database. Hazard ratios (HR) and 95% confidence intervals (CI) were calculated in this study.

The nomogram was used to predict cancer prognosis individually by incorporating clinical characteristics and risk scores of the patients. The calibration curves were utilized to visualize the deviation of predicted probabilities from what actually happened. The concordance index (C-index) was applied to measure the predictive accuracy of the nomogram. Time-dependent Receiver operating characteristic (ROC) curve was generated using R package survival ROC.

Tumor Immune Estimation Resource (TIMER 2.0) database (http://timer.cistrome.org/) was used to analyze the correlation between MD2 expression and the infiltration of six types of immune cells (B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages and dendritic cells) in LGG and GBM (39). In addition, CIBERSORT (http://cibersort.stanford.edu) was also applied to analyze the relationship between MD2 expression and 22 types of human immune cell subpopulations based on CGGA and TCGA datasets (40). The correlation between MD2 methylation and immune cell infiltration was analyzed by R package EpiDISH at cg13213009 and cg23732024 CpG sites. Estimation of stromal and immune cells in malignant tumor tissues using expression data (ESTIMATE) (https://bioinformatics.mdanderson.org/estimate) was employed to calculate the degree of immune cell infiltration.

KEGG and GO were applied to assess MD2 associated potential functions in gliomas based on TCGA database with R package ClusterProfiler (41). The correlation between MD2 expression level and immunomodulators was evaluated by TISIDB (http://cis.hku.hk/TISIDB) database in LGG and GBM respectively (42). The 50 interacting proteins with MD2 were collected from STRING (STRING: string-db.org) and top-100 MD2-related genes were obtained from GEPIA2 (gepia2.cancer-pku.cn) (43). The protein-protein interaction (PPI) network was visualized by Cytoscape software (44). Pearson’s correlation analysis was conducted between MD2 and the coincide of interaction and related genes using the GEPIA2 in LGG and GBM respectively.

The glioma cell lines U87 and A172 were purchased from the Chinese Academy of Sciences Cell Bank (Shanghai, China) and cultured in Dulbecco’s Modified Eagle Medium (DMEM) containing 10% fetal bovine serum (FBS), penicillin (100 U/mL) and streptomycin (100 μg/mL) and cultured and humidified incubator which maintained at 5% CO₂ and 37°C.

The glioma cell lines were transfected with 100 nM siRNA targeting MD2 (sense, 5’-GAAUCUUCCAAAGCGCAAATT-3’) or a non-coding scramble negative control siRNA (sense, 5’- TTCTCCGAACGTGTCACGTTT-3’) using 3 μL RNAi MAX reagent (Invitrogen, USA) in the opti-MEM medium. After 6 h incubation, the media was changed to normal DMEM and then cultured for 48 h.

Cells were lysed in NP-40 lysis buffer (150 mM NaCl, 100 mM NaF, 50 mM Tris-HCl (pH 7.6), and 0.5% NP-40) with protease inhibitor PMSF (1:100) at 4°C for 15 min. The protein concentration was measured using the BCA protein assay kit (Beyotime, Shanghai, China). The protein samples were subjected to 12% SDS-PAGE and then transferred to PVDF membrane and blocked in TBST with 5% no-fat milk at room temperature for 2 h. The membranes were incubated with primary antibodies at 4°C overnight. Then, the membranes were washed in TBST buffer and incubated with species-match HRP-linked secondary antibody at room temperature for 2 h. Afterwards, membranes were washed three times in TBST buffer, developed using the enhanced chemiluminescence (ECL) reagent (Fdbio science, Hangzhou, China), and captured by Tanon-5200Multi Imaging System (Tanon, Shanghai, China). Antibodies used in this study were as follows: rabbit polyclonal anti-MD2 (1178-1-AP; Proteintech, Chicago, USA), mouse monoclonal anti-β-Tubulin (2128; Cell Signaling Technology, Boston, USA), anti-mouse IgG, HRP-linked antibody (7076S; Cell Signaling Technology, Boston, USA) and anti-rabbit IgG, HRP-linked antibody (7074S; Cell Signaling Technology, Boston, USA).

Total RNA was extracted from the cells using the TRIzol reagent (Accurate Biology Co. Ltd, Hunan, China) according to the manufacturer’s instruction. Then, cDNA was synthesized by using reverse transcriptase (Vazyme, Nanjing, China). qRT-PCR was carried out using SYBR green supermix (Vazyme, Nanjing, China) with gene-specific primers. β-Actin was used as an internal standard for normalization. The sequences of primers for qRT-PCR are provided in the Table 1 . Each assay was performed in triplicate and the data were analyzed with the 2^-△△Ct.

The R software (Version 4.1.0), Graphad Prism 8 software (Version 8.0.2) and Adobe Illustrator software (Version 24.0.2) were used to perform statistical analysis and generate figures. Difference analysis between two groups was analyzed using Student’s t-test, and p < 0.05 was considered statistically significant. Continuous variable fitting a normal distribution was described as the mean with standard deviation.

The expression level of MD2 in glioma tissues and normal brain tissues was analyzed using the data from CGGA databases. The analysis indicated that MD2 was significantly overexpressed in glioma patients in comparison to normal brain tissues ( Figure 1A ). Further analyses of other two databases of TCGA and Rembrandt obtained similar results ( Figures 1B, C ). Next, to analyze the association between the expression level of MD2 and overall survival of glioma patients, Kaplan-Meier plotter analysis was performed with the datasets from CGGA, TCGA and Rembrandt. According to the median expression of MD2 in each dataset, the expression level of MD2 in glioma patients was separated into two groups with high and low expression. The Kaplan-Meier curves showed that high expression level of MD2 was remarkably related to the poor overall survival of glioma patients in CGGA ( Figure 1D , p < 0.0001), TCGA ( Figure 1E , p < 0.0001) and Rembrandt ( Figure 1E , p < 0.0001), respectively. These results indicated that MD2 could function as an oncogene, and its high expression may portend a worse prognosis in gliomas.

To elucidate potential roles of MD2 in the malignant progression of gliomas, we analyzed its expression levels in different grades of glioma in the datasets of CGGA and TCGA. Although there was no significant difference of MD2 expression between grade II and III in the dataset of CGGA, MD2 expression level was significantly increased along with the progression of gliomas from grade II to grade IV in both datasets ( Figure 2A ). Since IDH1 mutation is recognized as a principal driver in low grade gliomas, with an incidence of more than 70% (45, 46), we therefore examined the relationship between MD2 expression and the status of IDH1. In both databases from CGGA and TCGA, patients with higher MD2 expression level were synchronized with wild-type IDH1, whereas most of those with lower MD2 expression was associated with IDH1 mutation ( Figure 2B ). In parallel, 1p/19q codeletion is an important clinicopathologic characteristic for gliomas progression, and codeleted patients usually survive longer than non-codeleted patients (46, 47). Thus, we assessed the potential clinical association between MD2 expression and the status of 1p/19q in both databases of CGGA and TCGA. The analyses indicated that MD2 expression was significantly upregulated in 1p/19q non-codeleted group compared to the patients with 1p/19q codeleted ( Figure 2C ). We further analyzed the expression of MD2 in glioma patients with different ages in the databases of CGGA and TCGA. According to the median ages of patients, the glioma patients were separated into high-age (> 43 or > 46 years) and low-age (≤ 42 or ≤ 46 years) groups. The results indicated that the expression level of MD2 was significantly lower in the low-age group than that in the high-age group ( Figure 2D ). In addition, we found that the expression level of MD2 was increased after chemotherapy and radiotherapy, suggesting that the expression of MD2 may be related to therapeutic resistance ( Figures 2E, F ). These results indicated that high expression of MD2 was correlated with faster progression of gliomas.

To investigate the causes of abnormal expression of MD2 in gliomas, we detected MD2 expression level and its DNA methylation status. As shown in Figure 3A , we observed an obviously negative correlation between MD2 expression and its DNA methylation at two CpG sites including cg13213009 (R = -0.4903, p < 0.0001) and cg23732024 (R = -0.4499, p < 0.0001), while the methylation at cg17503786 did no correlate with MD2 expression. Subsequently, we selected the sites of cg13213009 and cg23732024 CpG to further examine the prognostic values of MD2 methylation in glioma patients. Consequently, MD2 methylation levels at both CpG sites significantly decreased in accordance with the progression of gliomas from grade II to grade IV ( Figure 3B ). Next, we established the relationship between MD2 methylation level and the status of IDH1. The results showed that higher methylation levels at both CpG sites tend to be associated with IDH1 mutation ( Figure 3C ). Moreover, Kaplan-Meier plots indicated that lower DNA methylation of MD2 gene correlates with shorter overall survival of glioma patients (p < 0.001) based on the analyses of CGGA and TCGA databases ( Figure 3D ). Collectively, these analyses demonstrated that MD2 gene methylation was also associated with the progression of glioma and the increased expression level of MD2 in glioma was induced by its reduction of DNA methylation.

To explore whether MD2 is an independent and significant factor for the prognosis of gliomas, we carried out univariate and multivariate Cox regression analyses using the data from CGGA and TCGA. The results indicated that MD2 expression (univariate hazard ratio (HR): 0.6, p = 1.4e-07; multivariate HR: 0.8, p = 4.2e-02), 1p/19q codeletion (univariate HR:0.3, p = 1.2e-12; multivariate HR: 0.4, p = 1.4e-05) and IDH1 mutation (univariate HR:0.3, p < 2e-16; multivariate HR: 0.5, p = 5.2e-07) could serve as independent protectable variances for gliomas, and WHO grade (univariate HR:2.8, p < 2.6e-16; multivariate HR: 2.8, p = 2.1e-34) and age (univariate HR:1.7, p = 1.2e-06; multivariate HR: 1.4, p = 8.4e-03) could be risk factors ( Figure 4A ). Similar results were also obtained by using the data from TCGA ( Figure 4B ). Next, we constructed the nomograms with these independent prognosis factors (Age, WHO grade, IDH1 status, 1p/19q status and MD2 expression) to predict 1-, 3- and 5-year survival probability of each glioma patient ( Figure 5A ). The calibration plot for the probability of survival revelated that the nomogram-predicated survival probability was very close to the ideal reference line in the databases of CGGA and TCGA, and the C-index were 0.75 in CGGA dataset and 0.86 in TCGA dataset, respectively ( Figure 5B ). In addition, the areas under ROC curve (AUC) of MD2 expression is 0.643 in 1-year survival, 0.675 in 3-year survival and 0.668 in 5-year survival in CGGA database and is 0.781 in 1-year survival, 0.762 in 3-year survival and 0.687 in 5-year survival in TCGA database ( Figure 5C ). Since WHO grade, IDH mutant status and 1p/19q codeletion are important clinicopathologic characteristics for glioma progression, we further performed ROC analysis combining MD2 expression with these parameters. As a result, the AUC of 1-, 3-, and 5-year survival rates in CGGA database were 0.643, 0.675 and 0.668, respectively. Similar results were also obtained in TCGA database ( Figure 5C ). These results suggested that MD2 is an independent factor for the prognosis of glioma.

To dissect the biological functions of MD2 in gliomas, we performed GO and KEGG enrichment analyses based on MD2 expression. The GO functional analyses indicated that a variety of functions are associated with MD2, among which the major functions are linked to T-cell costimulation, innate immune responses and inflammatory responses ( Figure 6A ). Meanwhile, KEGG pathway analyses showed that MD2-related pathways involve neuroactive ligand-receptor interaction, phagosome, infection and leishmaniasis ( Figure 6B ). To clarify these results, we firstly examined MD2-related immune-inhibitors and immune-stimulators using TISIDB database. We found that 8 immunoinhibitors (CD96, CSF1R, HAVCR2, IL10, IL10RB, LGALS9, PDCD1LG2, and TGFBR1, R > 0.5) and 6 immunostimulators (CD28, CD40, CD48, CD86, IL2RA, TMEM173, R > 0.5) were significantly associated with MD2 in gliomas ( Figure 6C ). Then, we detected the association between MD2 expression and nine immune checkpoints (PD1, PDL1, PDL2, LAG3, CTLA4, TIGIT, IDO1, CD276, CD47), which are promising immunotherapeutic targets for gliomas (48, 49). The analysis revealed that MD2 is positively associated with PDL1, PDL2 and CD276 ( Figure 6D ). To further analyze the potential roles of MD2 in glioma progression, MD2-binding proteins and genes correlated with MD2 expression were identified using the databases of STRING and GEPIA2. Fifty MD2-binding proteins and top 100 genes correlated with MD2 were obtained and the PPI network of those genes was mapped as shown in Figure 6E . The Venn diagram indicated that four common members, including CD14, LY86, TLR1 and TLR4, occur in both groups ( Figure 6F ), suggesting their regulatory roles in mediating the innate immune responses and inflammatory responses. The subsequent correlation analyses also indicated MD2 was strongly correlated with CD14, LY86, TLR1 and TLR4 in gliomas ( Figure 6G ). Through these functional and pathway enrichments, we found that MD2 strongly correlates with immunological responses in gliomas.

Given the involvement of MD2 in immunomodulatory signaling pathways in gliomas, we explored the association between MD2 expression and immune cell infiltration. The TIMER algorithm was used to examine the correlation between MD2 expression and six types of immune cell infiltration in LGG and GBM, respectively. MD2 showed a significantly positive correlation with the infiltrations of macrophages, neutrophils and NK cells in both LGG and GBM ( Figure 7A ). Meanwhile, the immune scores were calculated by ESTIMATE database, showing that immune scores were positively related with MD2 expression (p < 0.0001) ( Figure 7B ). Additionally, the correlation between MD2 and 22 types of infiltrating immune cells was calculated by CIBERSORT algorithm using the data from CGGA and TCGA databases. Different from the analysis results by TIMER, in addition to macrophages and neutrophils, three other immune infiltrating populations, including B cells, CD4⁺ T cells and CD8⁺ T cells, were also significantly and positively correlated with MD2 expression ( Figure 7C and Table 2 ). Furthermore, we found that DNA methylation of MD2 gene was negatively correlated with the infiltration of neutrophils, which is consistent with our analysis ( Figure 7D ). The correlation between MD2 expression and classical phenotypes of macrophages and neutrophils were analyzed in the databases of CGGA and TCGA, and we discovered that MD2 possesses an exceptionally positive correlation with M0 and M2 markers of TAMs (M2-type macrophages promote tumor progression) instead of M1 marker (M1-type macrophages inhibit tumor progression) ( Figure 7E ). Similarly, MD2 also showed a positive correlation with N2 phenotype marker of TANs ( Figure 7F ). The correlations between MD2 expression and the expression of other immune cell-specific markers were further confirmed by analyzing the database of TIMER as shown in Table 3 .

To ensure the correlation between MD2 expression and the infiltrations of macrophages and neutrophils, MD2 was knocked-down using siRNA in U87 and A172 glioma cells, resulting in remarkable decrease of the expression of MD2 at mRNA and protein levels ( Figure 7G ). Subsequently, we detected the changes of mRNA expression of the cytokines secreted by glioma cells. Both M2-type TAMs-related polarization factors (CSF-1, CCL-2, IL-10 and TGF-β) and N2-type TANs-related polarization factors (CXCL-2, CXCL-5, G-CSF and GM-CSF) were markedly decreased after the knockdown of MD2 in A172 and U87 cells ( Figure 7G ). At the same time, mRNA level of TLR4 was also dramatically reduced by silencing MD2 in A172 and U87 cells ( Figure 7G ). Collectively, these data revealed that MD2 expression is significantly associated with macrophage and neutrophil infiltration to promote their polarization toward M2-TAMs or N2-type TANs.