The overview of our study design is shown in Supplementary Figure S1. The LGGs gene expression data and matched clinical annotations were retrieved from TCGA and Chinese Glioma Genome Atlas (CGGA) databases. Only patients with survival time greater than 30 days were involved in this study. In total, four cohorts were enrolled: TCGA-LGG (481 samples), CGGA-mRNAseq_693 (420 samples), CGGA-mRNAseq_325 (170 samples), and CGGA-mRNA-array_301 (159 samples). The TCGA-LGG cohort data by log2(x+1) transformed were downloaded from the UCSC Xena website (https://xena.ucsc.edu/), and the other three cohorts were downloaded from the CGGA official website (http://www.cgga.org.cn/index.jsp). At the same time, the mRNA sequencing (mRNA-seq) data in CGGA was also log2(x+1) normalized.

We collected the genes of the main glycosyltransferase families involved in human species-specific glycosylation mechanisms (Mohamed Abd-El-Halim et al., 2021), including 221 genes retrieved from the GlycoGene Database (GGDB, http://riodb.ibase.aist.go.jp/rcmg/ggdb) and 213 genes extracted from the Hugo Gene Nomenclature Committee (HGNC, https://www.genenames.org/data/genegroup/#!/group/424). Eventually, a total of 143 overlap genes were selected in our analysis. Consensus clustering (Zhang et al., 2020) was performed to identify distinct glycosylation subtypes according to the expression of these 143 genes in the TCGA cohort. The ConsensuClusterPlus package was used and 1000 times repetitions were conducted to ensure the stability of clustering results. In addition, we applied the T-distributed stochastic neighbor embedding (t-SNE) to verify the reliability of subtype assignments based on the expression of the aforementioned genes (Bai et al., 2021).

To explore the variations in biological processes between glycosylation modification patterns, we selected the GSVA enrichment analysis to accomplish this purpose by using the “GSVA” package. GSVA is usually applied to estimate the variation of pathway and biological process activity in a nonparametric and unsupervised manner (Hänzelmann et al., 2013). The gene sets of “c2.cp.kegg.v7.4.symbols”, GOBP, and metabolism-associated pathway signatures derived from KEGG and REACTOME were both downloaded from the molecular signature database (MsigDB, http://www.gsea-msigdb.org/gsea/msigdb/search.jsp). Adjusted p < 0.05 was considered statistically significant.

The single sample gene set enrichment analysis (ssGSEA) was applied to estimate the enrichment scores of 13 immune functions and 16 immune cells. The scores were normalized to specific distribution from 0 to 1 and could quantify the relative abundance of estimated items. The marker gene sets were curated from a current study (Bindea et al., 2013), and the “GSVA” package was introduced to conduct the aforementioned analysis. Moreover, to comprehensively estimate the infiltration level of immune cells and stromal cells in the TME, Immune Cell Abundance Identifier (ImmuCellAI) (Miao et al., 2020), Tumor Immune Estimation Resource (TIMER) (Yu L. et al., 2021), and microenvironment cell population counter (MCP-counter) (Ye et al., 2021) algorithms were also utilized in our study. Among them, the former is calculated in the web tool ImmuCellAI (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/), while the latter two algorithms are implemented on the Sangerbox platform (http://vip.sangerbox.com/home.html).

To identify glycosylation-related genes (GRGs), we performed differential expression analysis between glycosylation distinct subtypes. The “limma” package was used to determine differentially expressed genes (DEGs), and the significance filtering criteria were set as | log2 (fold change) | > 1 as well as adjusted p value less than 0.0001. In addition, to ensure the reliability of GRGs, weighted gene co-expression network analysis (WGCNA) was also introduced to establish co-expressed gene modules strongly related to glycosylation. Based on the expression of 8619 prognostic genes (univariate Cox analysis p < 0.01), a gene co-expression network was established by using the “WGCNA” R package (Langfelder and Horvath, 2008). The Pearson correlation was calculated according to the expression of these filtered genes, and an optimal soft-thresholding power β = 7 was selected to build an unsigned weighted adjacency matrix followed by a topological overlap matrix (TOM) conversion. Finally, the average linkage hierarchical was clustered with the parameter height = 0.25, and the criterion for gene module identification was set as a cutting height of 0.9 along with the module genes a minimum number of 20. The modules of | correlation coefficient | > 0.5 were considered as closely related to glycosylation and selected for a subsequent analysis. The intersections of DEGs and module genes screened by WGCNA were defined as GRGs.

To quantify the glycosylation modification patterns of individual patients, we developed a glycosylation scoring scheme to assess the glycosylation modification pattern of individuals with LGGs, which was glycosylation gene signature and the glycosylation score. Primarily, we divided the GRGs into two groups based on their hazard ratio (HR) of the univariate Cox regression model. Genes with HR > 1 and those with HR < 1 were defined as the glycosylation gene signatures A and B, respectively. Next, we conducted GSVA for dimension reduction of the signatures A and B, the GSVA value of which were separately termed as value A and value B. Finally, both value A and value B were selected to calculate the glycosylation score as follows: glycosylation score = value A- value B. Patients with LGGs were classified into high and low glycosylation score groups based on the median cutoff value, and scores between the glycosylation subtypes or clusters were assessed by the Wilcoxon test.

Moreover, the prognostic value of the glycosylation scores was validated in three cohorts of the CGGA database by using the same signatures and the cutoff value.

To explore the association between the glycosylation score and some related biological pathways. The gene sets related to stromal activation were extracted from the supplementary material of a current study (Mariathasan et al., 2018), including epithelial–mesenchymal transition (EMT) signatures (EMT1, EMT2, and EMT3) and pan-fibroblast TGFb response signature (Pan-FTBRS). Moreover, immune checkpoint- and immune activation–related genes were retrieved from the publication of Zhang et al. (2020). In addition, immunosuppressive factors and immunosuppressive cell recruitment factors were obtained from another study (Su et al., 2019). The differences in biological process between distinct glycosylation score groups were assessed via the Wilcoxon test.

Differentially expressed genes (DEGs) between the high and low glycosylation score groups were set as adjusted p-value < 0.001 and | log2 (fold change) | > 1, which were determined by using the “limma” package. The DEGs were visualized into volcano plots and submitted to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses by the “clusterProfiler” package. Furthermore, both of the top 500 up- and downregulated genes were uploaded into the Connectivity Map 02 (CMap, https://portals.broadinstitute.org/cmap/) database, and a CMap mode-of-action (MoA) analysis was applied to discover possible small molecular drugs for LGGs and mechanisms of action.

We evaluate the predictive capacity of glycosylation score in responding to chemotherapeutic agents and immunotherapy. The 50% inhibiting concentration (IC50) value of 138 drugs of each patient was quantified using the pRRophetic algorithm (Tan et al., 2021), and the differences between high and low glycosylation score subgroups compared by using the Wilcoxon test. In addition, the response to immunotherapy was predicted by using immunophenoscore (IPS) (Gui et al., 2021; Wu et al., 2021), which was obtained without bias based on determining components of immunogenicity including effector cells, immunosuppressor cells, major histocompatibility complex (MHC) molecules, and immune modulators. The statistical difference of the IPS value was tested by the Wilcoxon test.

R software (version 4.1.0) was applied to conduct statistical analysis. We adopted the Kaplan–Meier (K-M) analysis and log-rank tests to perform the survival analysis with the “survival” R package. The prognostic value of glycosylation score was evaluated by a multivariate Cox regression model and a time-dependent receiver operating characteristic (ROC) analysis using R packages “survival’ and “timeROC”. Comparisons between glycosylation subgroups, clusters, and score groups were presented via the Wilcoxon test.

Based on the expression of 143 genes involved in glycosyltransferase families, 481 patients from the TCGA database were obviously divided into two groups when k = 2 (Figure 1A) by the consensus clustering analysis, which was termed as glycosylation subtype A and subtype B, respectively. The t-SNE algorithm confirmed that glycosylation subtype assignments are reliable and can be completely distinguished (Figure 1B). The survival analysis revealed that patients of subtype A showed a significantly shortened survival time compared with patients of subtype B (Figure 1C, log-rank test).

To explore the biological molecular changes between distinct glycosylation subtypes, the KEGG gene set was used to performed GSVA enrichment analysis. As shown in Figure 1D, subtype A presented enrichment pathways prominently related to carcinogenic activation and immune regulation pathways such as nucleotide excision repair, homologous recombination, antigen processing and presentation, and NOD-like receptor signaling pathway, while subtype B is markedly enriched in metabolism-related processes, including steroid biosynthesis, oxidative phosphorylation, citrate cycle TCA cycle, and so on. Next, we further compare the metabolism and immune microenvironment features affected by glycosylation modification patterns between the subtypes. The results of ssGSEA showed a significant up-regulation in immune infiltrating cells (Figure 1E) and immune function (Figure 1F) of subtype A. Furthermore, the GSVA results revealed significant differences in metabolic processes (Figure 2A), with subtype B mainly active in lipid metabolism (Figure 2B) and other metabolism (Figure 2C), while subtype A was mainly active in carbohydrates metabolism (Figure 2D) and amino acid metabolism (Figure 2E).

To further validate the glycosylation modification patterns and the potential biological behavior, GRGs identified by differentially expressed analysis and WGCNA were further inputted into the subsequent investigation. Primarily, 1945 DEGs between glycosylation subtypes including 788 up-regulated genes and 1157 down-regulated genes were visualized in the volcano plot (Figure 3A), and the top 100 ones are displayed in Figure 3B. Then 18 gene modules (Figure 3C) were recognized by WGCNA, and six modules (light yellow, brown, grey60, dark orange, dark grey, and tan) were closely related to glycosylation (Figure 3D). Eventually, 971 interactions (Figure 3E), considered GRGs (Supplementary Table S1), were uploaded into the consensus clustering analysis. Consistent with the subtype grouping of glycosylation modification patterns, the TCGA-LGGs cohort was divided into two distinct clusters, which were defined as glycosylation clusters A and B (Figure 3F). Also, the t-SNE algorithm data confirmed the rationality of the cluster assignments (Figure 3G). As expected, we found that patients in cluster A exhibited significantly shortened the survival time (Figure 3H; log-rank test). This demonstrated that two distinct glycosylation modification patterns did exist and work in LGGs. We further probed into the metabolism and immune microenvironment characteristics of the glycosylation clusters. In short, glycosylation cluster A showed higher scores of immune infiltrating cells (Figure 4A) and immune function (Figure 4B), except for CD8⁺ T cells, mast cells, and natural killer (NK) cells. Moreover, metabolic characteristics also showed significant differences (Figures 4B–E) between glycosylation clusters. Excitedly, these results were much similar to that of between glycosylation subtypes. Collectively, the concordance between prognostic, metabolic, and immune microenvironment characteristics in glycosylation clusters and subtypes confirmed again that glycosylation modification played a vital role in metabolic reprogramming and immune regulation, and this classification was steady and reliable.

Although our findings revealed the role of glycosylation modification in prognosis, metabolic reprogramming, and immune modulation, these analyses could not accurately quantify the glycosylation modification patterns of individual tumors. Therefore, to enable the quantification of glycosylation patterns in individual LGGs patients, we developed a glycosylation scoring system termed as the glycosylation score. Based on the median glycosylation score of 0.009856, patients were classified into the high and low glycosylation score subgroups. Figure 5A showed the distribution of clinical traits including age, gender, grade, overall survival, and survival status between high and low glycosylation score subgroups, while Figure 5B further displayed the grade distribution, glycosylation subtypes, and clusters assignment of patients from distinct glycosylation score subgroups, which are summarized in Supplementary Table S2. The Wilcoxon test demonstrated that patients in subtype A had a higher glycosylation score compared to subtype B (Figure 5C); meanwhile, cluster A showed a higher glycosylation score than that of cluster B (Figure 5D). This indicated that a high glycosylation score could be closely related to poor prognosis, whereas a low glycosylation score could be closely related to a favorable prognosis. Of course, this hypothesis was confirmed by a subsequent survival analysis (Figure 5E) based on glycosylation scores. We analyzed the prognostic values of glycosylation scores in depth. The multivariate cox regression analysis demonstrated that the glycosylation score was an independent risk factor with its HR = 2.010 for LGGs (Figure 5F). Meanwhile, the time-dependent area under the curve (AUC) suggested that the glycosylation score had a robust and reliable value in predicting the prognosis for LGGs in the TCGA cohort (Figure 5G).

We further verified the predictive effectiveness of glycosylation score in three independent external cohorts according to the same glycosylation gene signatures and cutoff score. Consistent with findings in the TCGA cohort, patients with a high glycosylation score showed a significantly shortened survival than those with a low glycosylation score in validation cohorts including CGGA-mRNAseq_693 (Supplementary Figure S2A), CGGA-mRNAseq_325 (Supplementary Figure S2B), and CGGA-mRNA-array_301 (Supplementary Figure S2C). Similarly, superior time-dependent AUCs were observed in those cohorts (Supplementary Figures S2D–F).

The correlations between the glycosylation score and metabolism as well as immune microenvironment features were further explored in this study. Consistent with findings observed in distinct glycosylation patterns, the glycosylation score was strongly linked to metabolism activity (Supplementary Figure S3A) and immune microenvironment infiltrations (Figure 6). The results showed that the high glycosylation score subgroup was associated with enhanced activity in most amino acid metabolism (Supplementary Figure S3B) and carbohydrates metabolism (Supplementary Figure S3C); on the contrary, the low glycosylation score subgroup was associated with enhanced activity in most lipid metabolism (Supplementary Figure S3D) and other metabolism (Supplementary Figure S3E). Meanwhile, the distribution of infiltrating cells estimated by ImmuCellAI, TIMER, MCP-counter, and ssGSEA algorithms between the high and low glycosylation subgroups were also investigated. Our results (Figure 6) suggested that the high glycosylation score subgroup had a significantly higher infiltrating level than the low glycosylation score subgroup. Most of the infiltrating immune and stromal cells increased in the high glycosylation score subgroup, whereas T helper 1 (Th1) and CD8⁺ naïve cells decreased in this subgroup. Furthermore, there were higher scores of 13 immune functions in the high glycosylation score subgroup (Figure 7A), compared to the low glycosylation score subgroup. Taken together, the glycosylation score presented a close association with metabolism and immune microenvironment features of LGGs.

To explore potential mechanisms correlated with different glycosylation score subgroups in depth, we compared the expression of chemokine and cytokine between subgroups. These cytokine and chemokine retrieved from published literature, of which, includes the transcripts correlated with immune activation (Figure 7B), immunosuppressive cell recruitment factors (Figure 7C), immunosuppressive factors (Figure 7D), and immune checkpoints (Figure 7E). We observed that most of these genes were expressed higher in the high glycosylation score subgroups, revealing the coexistence of a chronic inflammation state and a suppressed immune microenvironment in this subgroup. To better depict the role of glycosylation scores, we evaluated other related biological processes, including immune activation, stromal activation, and cancer promotion, in patients with LGGs. Consistent with the aforementioned observations, patients in the high glycosylation score subgroup exhibited a much higher enrichment score of immune activation–related biological processes (Figure 7F). On the other hand, higher stromal activation and cancer promotion scores were also enriched in this subgroup. Collectively, the glycosylation score did play a non-negligible role in affecting LGG progression.

To explore potential small therapeutic drugs for patients with high glycosylation scores, 2090 DEGs, including 923 up-regulated genes and 1167 down-regulated genes, between the high and low glycosylation score subgroups were identified (Figure 8A). The GO analysis results illustrated that the DEGs are primarily involved in modulation of chemical synaptic transmission, regulation of trans-synaptic signaling, synaptic membrane, ion channel complex, gated channel activity, and ion channel activity (Figure 8B). Meantime, these DEGs were significantly enriched in the cholinergic synapse, GABAergic synapse, glutamatergic synapse, neuroactive ligand–receptor interaction, serotonergic synapse, and synaptic vesicle cycle signaling pathways (Figure 8C). Based on the CMap MoA analysis, a total of 14 potential small therapeutic drugs (such as spiradoline, propofol, and dextromethorphan) and 13 drug mechanisms (such as opioid receptor agonist, GABA receptor agonist, and glutamate receptor antagonist) were identified (Figure 8D).

To explore the potential role of glycosylation score in the drug sensitivity prediction, we compared the IC50 value of 138 drugs in TCGA-LGGs patients by using the pRRophetic algorithm; of which, the results showed that there were 29 drugs (such as gefitinib, nilotinib, and axitinib) in which the estimated IC50 value were significantly lower in the low glycosylation score subgroup compared to the high glycosylation score subgroup (Supplementary Figure S4A), indicating that patients in the low glycosylation score subgroup might respond better to these chemotherapy drugs. Similarly, patients in the low glycosylation score subgroup had a higher IPS score, suggesting that patients in this subgroup might be more sensitive to immunotherapy (Supplementary Figure S4B).