All the information of patients with grade II–III glioma was obtained from The Cancer Genome Atlas database (Center for Cancer Genomics at the National Cancer Institute; http://portal.gdc.cancer.gov/). The Perl language (http://www.perl.org/) was used to merge RNA-sequencing results into matrix files. The data of 479 patients were obtained and the patients were randomly divided into two groups using the R (version 4.2.3; http://www.R-project.org/) package ‘caret’ (20). The first group (n=240) was used for training and the other group (n=239) was used for testing. All data regarding the clinical characteristics of all patients are shown in Table SI.

A total of 243 glycolysis-related gene profiles (21) (Table SII) were obtained from the Molecular Signatures Database (http://www.gsea-msigdb.org) labeled ‘Hallmark glycolysis’ and ‘Kyoto Encyclopedia of Genes and Genomes (KEGG) glycolysis gluconeogenesis’.

The lncRNA annotation files were obtained from the Generic Code Database (https://www.gencodegenes.org/human/). Subsequently, lncRNA information related to glycolysis-related gene expression was obtained. The selection criteria were set as |R|>0.5 and P<0.001.

The training set was used to construct the prognostic model. As previously reported by Sun et al (22), a gene signature with high performance and stability was constructed using a combination of machine learning algorithms, including the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm, and multivariate Cox analyses. Univariate Cox analysis was applied to determine which glycolysis-related lncRNAs were associated with overall survival (OS) among patients with grade II–III glioma (P<0.05). Subsequently, LASSO regression analysis was conducted to analyze the glycolysis-related lncRNAs using the package ‘glmnet’ in R (23). Upon reaching the minimum partial likelihood deviance, the optimal lambda was identified. Finally, a risk characteristics and prognosis model based on seven glycolysis-related lncRNAs was constructed as follows: Risk score=∑k=1ncoef(incRNA)*expr(incRNAk)

Coef (lncRNA) refers to the coefficient associated with survival. Expr is the expression of lncRNAs.

Kaplan-Meier (K-M) analysis was applied to compare OS between two groups using the R package ‘survival’ (24) (version 3.5–0). The package ‘timeROC’(25) (version 0.4) was used to verify the predictive ability of the established risk model. The concordance index (C-index) was estimated to evaluate the discrimination of the model. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis were used to visually distinguish two groups.

Independent prognostic factor analysis was conducted in R with the package ‘survival’. The R package ‘rms’ (26) was further applied to construct a nomogram and the predictive ability of the nomogram was evaluated.

To quantitatively analyze tumor tissue transcriptomic data, the CIBERSORT algorithm was used, which led to the estimation of the absolute number of immune and stromal cells present in the tumor tissues. Furthermore, the single sample Gene Set Enrichment Analysis (ssGSEA) algorithm was used to assess the difference between high- and low-risk groups in terms of the ability to resist tumor infiltration (27). The ESTIMATE algorithm was used to determine the stromal score, immune score and ESTIMATE score of grade II–III glioma samples (28). The package ‘maftools’ was used to analyze the TMB (29).

Drug information was obtained from the Genomics of Cancer Drug Sensitivity (GDSC) website (30) (http://www.cancerrxgene.org), and the R package ‘pRRophetic’ (31) was used to predict IC₅₀ values to explore possible clinical adjuvant drugs for the treatment of grade II–III glioma (https://www.cancerrxgene.org). To determine if the risk model was associated with immunotherapy, the expression levels of key genes in the two groups of immune checkpoints were compared.

Differentially expressed genes (DEGs) were identified using the package ‘limma’ (32). The ‘clusterProfiler’ package was used to study the related functions and pathways of DEGs bias (33). GSEA was used to compare potential pathways between the two groups. The R package ‘ggalluvival’ (34) (version 0.12.5) was used to visualize the association between lncRNAs and risk factors (protective/risk) in a Sankey diagram.

Cytoscape software (version 3.9.0; http://cytoscape.org) was used to visualize the lncRNA-mRNA co-expression network.

The HA1800 human normal astrocyte cell line and the SHG-44 human glioma cell line were purchased from American Type Culture Collection (35). HA1800 and SHG-44 cells were cultured in DMEM (Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (Gibco; Thermo Fisher Scientific, Inc.) and antibiotics (100 U/ml penicillin and 100 µg/ml streptomycin). Cells were incubated at 37°C with 5% CO₂.

Colorectal neoplasia differentially expressed (CRNDE) small interfering RNA (siRNA), si-CRNDE-1 (5′-GGTGTTAAGTGTGATGCTTCC-3′) or siCRNDE-2 (5′-GGATGCTGTCAGCTAAGTTCA-3′), and negative control siRNA (5′-UUCUCCGAACGUGUCACGUTT-3′) were purchased from Shanghai GenePharma Co., Ltd (cat. no. A10001). CRNDE expression in SHG-44 cells was silenced by siRNA transfection using Lipofectamine^® RNAiMax reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions. The transfection compound was prepared at room temperature for 20 min. The concentration of siRNA used was 50 nM. After culturing in the incubator (37°C) for 48 h, the cells and cell cultures were harvested for subsequent experiments, and the transfection efficiency was evaluated by reverse transcription-quantitative PCR (RT-qPCR).

The cell density for RNA extraction was 1×10⁶. Total RNA was extracted using the TRIzol^® reagent (Thermo Fisher Scientific, Inc.), and cDNA was obtained by reverse transcription using the TransScript^® First-Strand cDNA Synthesis SuperMix (Shanghai Yeasen Biotechnology Co., Ltd.). The cDNA was amplified using the Hieff UNICON universal Blue qPCR SYBR Green Master Mix kit (YEASEN Corporation, China). RNA extraction, cDNA synthesis and qPCR were all performed according to the manufacturer's protocols. The reaction volume was 20 µl, including 10 µl of Universal Blue qPCR SYBR Green Master Mix, 7.6 µl of nucleic acid-free water, 0.2 µl of each primer, and 2 µl of cDNA product. The PCR cycling conditions were as follows: 95°C for 2 min for 1 cycle, 95°C for 10 sec and 60°C for 30 sec for 40 cycles, and the melting curve was then determined. GAPDH was used as an internal reference and relative mRNA expression was calculated using the 2^−ΔΔCq method (36,37). The experiment was repeated three times. The primer sequences were as follows: GAPDH forward, 5′-AAGGTGAAGGTCGGAGTCAAC-3′ and reverse, 5′-GGGGTCATTGATGGCAACAAT-3′; and CRNDE forward, 5′-TGGATGCTGTCAGCTAAGTTCAC-3′ and reverse, 5′-TTCCAGTGGCATCCTCCTTATC-3′.

The effect of CRNDE expression on the proliferation of SHG-44 cells was detected using a CCK-8 assay (Dojindo Laboratories, Inc.). Briefly, after the cells were cultured in 96-well plates (6×10³ cells/well) for 0, 24, 48, 72 and 96 h, 20 µl CCK-8 cell proliferation agent was added to each well and incubated for 2 h at 37°C and the absorbance of the cell culture medium was measured at a wavelength of 450 nm.

After transfection, SHG-44 cells were cultured in 6-well plates (5×10³ cells/well) for 12 days. After fixation with 4% paraformaldehyde solution at room temperature for 20 min, the cells were stained with 0.1% crystal violet at room temperature for 15 min (Beyotime Institute of Biotechnology).

SHG-44 cells were seeded in 6-well plates (3×10⁵ cells/well) and incubated for 24 h. Prior to the assay, cells were subjected to serum deprivation to enhance glucose uptake. Subsequently, the cells were washed three times with PBS buffer. Following this, the cells were stimulated with 100 nM insulin for 30 min. A glucose uptake assay kit (Abcam) was used to quantify glucose uptake, with all results being normalized to the cell number.

All statistical analyses were conducted in the R and Perl software and the statistical results of samples were expressed as mean ± standard deviation. All assays were done at least three times independently. Independent Student's t-tests were applied to determine the difference between the two groups. Correlation analysis was performed using the Pearson correlation test. For the analysis of differences between K-M curves, the log-rank test was performed. If there is no special description for the above method, P<0.05 was considered to indicate a statistically significant difference.

Fig. 1 shows the process of the present study. A total of 14,056 lncRNAs and 243 glycolytic genes were screened for subsequent Pearson correlation analysis (21). The co-results among lncRNAs and mRNAs of Pearson correlation test are shown in Table SIII and 1,889 glycolytic lncRNAs were identified. Additionally, the consistent part of the expression between glycolytic genes and glycolytic lncRNAs was shown in the Sankey diagram, and their correlation was visually demonstrated (Fig. 2A and B).

Univariate Cox analysis was used to screen glycolytic lncRNAs. In the end, 656 prognostic lncRNAs were identified in the training set (Table SIV). LASSO regression analysis was performed on these survival-related lncRNAs to further screen key lncRNAs for model construction. Finally, nine lncRNAs were retained (Fig. 2C and D). Subsequently, multivariate Cox regression analysis was performed and seven glycolytic lncRNAs, including CRNDE, AC022034.1, Ras homolog family member Q (RHOQ-AS1), AL159169.2, AL133215.2, AC007098.1 and LINC02587, were identified (Fig. 2E; Table SV). Finally, the calculation formula was determined using the Cox regression model coefficient and lncRNA expression: CRNDE × 0.530754428201569 + AC022034.1 × 0.51597198629815 + RHOQ-AS1 × 4.0070016692408 + AL159169.2 ×-1.33589967405211 + AL133215.2 × 0.902923079138878 + AC007098.1 × 0.660037199843563 + LINC02587 × 0.29755586476417. AL159169.2, with a negative coefficient, tended to be a protective factor, and the other six lncRNAs tended to be risk factors, which was confirmed by the Sankey diagram (Fig. 3A). Additionally, the relative expression levels of the seven hub glycolytic lncRNAs were evaluated (Fig. 3B). Patients in the high-risk group had higher expression levels of destructive lncRNAs than patients in the low-risk group (Fig. 3A and B). Subsequently, the performance of the risk formula was verified using PCA and t-SNE algorithms. The results indicated that there was a difference between the two groups of patients with grade II–III glioma, which clearly showed that the situation of patients with grade II–III glioma in the two groups was different (Fig. 3C and D). Furthermore, the survival of the two groups of patients with grade II–III glioma in the training set were compared (Fig. 3E). The OS of different groups in the training set was compared by K-M analysis. The results demonstrated that patients with grade II–III glioma in the low-risk group had an improved OS (Fig. 3F; P<0.001). The receiver operating characteristic curve showed that the seven glycolytic lncRNAs had high accuracy in diagnostic efficacy, and the areas under the curve (AUCs) of 1-, 3- and 5-year prediction were 0.931, 0.921 and 0.866, respectively (Fig. 3G). The C-index of the risk model also showed that the risk model had good reliability (Fig. 3H).

The test set of patients with grade II–III glioma was subsequently used to evaluate the reliability of the model (n=239). K-M analysis indicated the survival of the two groups of patients with grade II–III glioma, and the results demonstrated that the OS of the low-risk group was improved compared with that of the high-risk group (Fig. 4A; P<0.001). Fig. 4B-D shows the distribution of the risk scores, survival time in the two groups and the expression of the seven lncRNAs in the test set. There was a notable difference between the two groups in the PCA and t-SNE analysis (Fig. 4E and F). In the test set, the predicted AUC values of the risk model for 1, 3 and 5 years were 0.846, 0.833 and 0.788, respectively (Fig. 4G). C-index analysis also illustrated that the risk model had good reliability (Fig. 4H). Finally, the entire set of patients with grade II–III glioma (n=479) was used to verify the prognostic ability of the model, indicating consistency with the aforementioned validation (Fig. 4I-P). It was evident that the established risk model was capable of reliably distinguishing patients with grade II–III glioma based on their risk profile.

To further estimate the distribution of the two groups, PCA was conducted using the whole gene expression profile, 243 glycolytic gene profiles, 1,889 glycolysis-related lncRNAs and the risk model (Fig. 5A-D). This risk model clearly indicated that patients with grade II–III glioma could be correctly classified by the metabolism of glycolysis between the two groups.

The next step was to verify the accuracy of the model based on seven lncRNAs associated with glycolysis for independent prediction of prognosis in grade II–III glioma. This was performed using univariate and multivariate Cox regression analysis. The results of univariate regression analysis [hazard ratio (HR), 1.003; 95% CI, 1.002–1.004; P<0.001] indicated that the model could accurately predict the prognosis of patients, and multivariate regression analysis (HR, 1.002; 95% CI, 1.001–1.003; P<0.001) also revealed similar results (Fig. 6A and B). In summary, the risk model was an independent prognostic factor for grade II–III glioma. A nomogram was constructed by combining multiple clinical factors with the constructed risk score to predict the survival rates of patients with grade II–III glioma (38) (Fig. 6C). Subsequently, the reliability of the nomogram was verified, and the results demonstrated that the actual observed value was similar to the predicted value (Fig. 6D). Additionally, the AUC values for the OS rate at each time point were 0.887, 0.884 and 0.825, respectively (Fig. 6E). Compared with the C-index of other clinical features, the nomogram and risk score had higher values (Fig. 6F).

In view of the extensive research on the TME in the evaluation and treatment of grade II–III glioma, a number of immune assessment methods were applied to determine the difference in the degree of immune infiltration between the two subgroups. First, some basic scores of patients with grade II–III glioma were analyzed using the ESTIMATE algorithm. The results of this indicated that with the increase of the risk level of patients with grade II–III glioma, the immune, stromal and ESTIMATE scores also increased (Fig. 7A-C). Using gene set variation analysis, it was revealed that high-risk grade II–III glioma was related to the change of some immune functions such as ‘cytolytic_activity’, ‘inflammation-promoting’ and ‘CCR’ (Fig. 7D). Subsequently, the categories and proportions of immune cells were analyzed using the CIBERSORT algorithm. The distribution differences and proportion changes of immune cells in the two groups were analyzed (Fig. 7E and F). Additionally, which immune cells were more common in high-risk populations was summarized (Fig. 7G. Finally, the immune cell infiltration of patients with grade II–III glioma in different risk groups was studied using the ssGSEA algorithm. The results indicated that the number of immune cell subsets in the high-risk group was higher (Fig. 7H). Additionally, some important immune functions were significantly upregulated in the high-risk group (Fig. 7I). In summary, it could be seen that patients with grade II–III glioma in the high- and low-risk groups had distinct immune infiltration differences, which may lead to differences in prognosis between the two groups.

TMB is a molecular marker that is used to quantify the mutations that tumor cells take as part of immunotherapy. It can evaluate whether immunotherapy is effective. The present study further explored the somatic mutation rate in the high- and low-risk groups. The results demonstrated that the mutation rate of the low-risk group [241 cases (97.57%) in 247 samples] appeared to be higher than that of the high-risk group [185 cases (85.65%) in 216 samples]. The top-ranked driver mutant genes are shown in Fig. 8A and B. Further quantitative analysis demonstrated that the TMB score was significantly increased in the high-risk group (Fig. 8C). Additionally, a positive correlation between risk score and TMB was demonstrated by Spearman correlation analysis. According to the median TMB score, grade II–III glioma samples were divided into two groups: High TMB group and low TMB group. According to the results of K-M analysis, the OS of the high TMB group was worse (Fig. 8E). When both risk score and TMB were considered in the K-M analysis, the risk score eliminated the better OS of low TMB. This suggested that the risk score could predict the prognosis of grade II–III glioma in a superior way than TMB (Fig. 8F). At the same time, the main immune cells associated with TMB were CD8 T cells and natural killer cells (Fig. 8G).

Considering the significant difference in prognosis between the two groups of grade II–III glioma samples, some potential drugs for precise treatment were selected. It was hypothesized that the response to drugs may be different between the two groups. The R package ‘pRRophetic’ was used to screen potential drugs based on the IC₅₀ value in the GDSC genomics database. The IC₅₀ estimates of A.770041, ABT.263, AP.24534, AS601245 and AZ628 were relatively high in low-risk samples.

This suggested that patients with grade II–III glioma in the high-risk group may benefit from the use of these compounds (Fig. 9A-J). The IC₅₀ values of axitinib and 5-aminoimidazole-4-carboxamide ribonucleotide were increased in the high-risk group, indicating that these two drugs may be used for the treatment of patients with grade II–III glioma in the low-risk group (Fig. 9F and G). The IC₅₀ of commonly used antitumor drugs in the two groups was determined. In the low-risk group, the IC₅₀ of cisplatin (Fig. 9H) and gemcitabine (Fig. 9I) was higher. The IC₅₀ of gefitinib was higher in the high-risk group (Fig. 9J), indicating that the risk score had high accuracy in predicting the sensitivity to antitumor drugs. Finally, it was explored whether the activation of immune checkpoint genes was consistent in the two groups. Consistent with our hypothesis, the activation of immune checkpoint genes in high-risk populations was generally relatively high (Fig. 9K and L).

In order to examine the mechanism that leads to significant differences between the two groups in a multi-level analysis, Gene Ontology (GO) and KEGG analyses were performed on 89 DEGs (|log₂ fold change|>2.0; P<0.05; Table SVI). According to GO analysis (Fig. 10A; Table SVII), DEGs were mainly enriched in tissue-related BPs. In terms of molecular functions, these DEGs were mainly concentrated in ‘extracellular matrix structural constituent’ and ‘glycosaminoglycan binding’. In addition, KEGG analysis showed that the DEGs were concentrated in ‘extracellular structure organization’ (Fig. 10B and C; Table SVII). Subsequently, through GSEA of KEGG, it was revealed that the pathways enriched in the high-risk (Fig. 10D; Table SVIII) and low-risk (Fig. 10E; Table SVIII) groups were different. Finally, a lncRNA-mRNA interaction network was constructed using Cytoscape software (Fig. 10F). The aforementioned results may provide some ideas regarding the mechanism of glycolysis-related lncRNAs in grade II–III glioma carcinogenesis.

Due to the extensive research on CRNDE, CRNDE was chosen for validation (39,40). The high CRNDE expression in grade II–III glioma was verified using cell lines. The experiments revealed high CRNDE expression in SHG-44 cells compared with HA1800 cells (Fig. 11A). In addition, the effects of CRNDE on cell functions were investigated. The function of CRNDE in SHG-44 cells was investigated via knockdown experiments and RT-qPCR confirmed that CRNDE was knocked down in SHG-44 cells (Fig. 11B). The results of the CCK-8 assay indicated that CRNDE could increase cell viability, which was the same as the results of database analysis (Fig. 11C). The colony formation assay also indicated that knockdown of CRNDE significantly inhibited the proliferation of SHG-44 cells (Fig. 11D and E). In addition, the potential regulatory role of CRNDE in glucose metabolism in grade II–III glioma was investigated using the ECAR method. The results demonstrated that knockdown of CRNDE reduced the ECAR and glycolysis (Fig. 11F and G). The experiments further confirmed the results of previous analyses.