A total of 700 grade II and III glioma samples were included in our study. The RNA-seq and clinical data were extracted from CGGA (http://www.cgga.org.cn/) and TCGA (https://portal.gdc.cancer.gov/) databases. In this study, TCGA-LGG dataset (n=522) was defined as the training cohort, whereas CGGA dataset (n=178) was the validation cohort. The characteristics of glioma samples in this study were summarized in Table 1 . Moreover, IMvigor210 dataset, a cohort of atezolizumab (anti-PD-L1 monoclonal antibody) for the treatment of urothelial carcinoma, was extracted to evaluate the predictive value of HRL signature for the efficacy of immunotherapy (16).

The 26 hypoxia-related genes were reported by previous studies (17, 18). A total of 14,488 and 13,895 lncRNAs were identified in the TCGA and CGGA datasets, respectively. Those lncRNAs whose expressions closely corelated with the expression of 26 hypoxia-related genes (|R|>0.5 and p<0.01) were identified as HRLs.

To evaluate the involvement of biological processes of each sample, gene set variation analysis (GSVA) was conducted to quantify the involvement of Gene Ontology pathways in each LGG sample (19). As for gene set enrichment analysis (GSEA), differentially expressed genes between two groups were identified. Those with false discovery rate (FDR) ≤0.05 were selected for GSEA analysis. Gene sets of hallmarks were obtained from the Molecular Signatures Database (MSigDB). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted by “clusterprofiler” R package. The nomogram and calibration curves were constructed and visualized using “rms” and “regplot” R packages. Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was used to predict the responses of glioma patients to immune checkpoint inhibitors (ICIs) (20).

Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) analysis was conducted to calculate the tumor purity of each sample by “estimate” R package (21). The infiltration of immune cells was estimated by single-sample Gene Set Enrichment Analysis (ssGSEA) and Tumor Immune Estimation Resource (TIMER), and Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithms (22–24).

The risk signature was constructed using the least absolute shrinkage and selection operator (LASSO) analysis. The risk score was calculated by the following algorithm:

For subgroup analysis, LGG patients were divided into different groups based on the following variables: grade (grades II or III), age (≤45 years old or >45 years old), and IDH status (mutant or wildtype).

The potential target miRNAs of LINC00941 and BASP1-AS1 were predicted using ENCORI online webtool (http://starbase.sysu.edu.cn/). Then, the potential target mRNAs of these miRNAs were predicted using microT-CDS (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index) and mirDIP online webtools (http://ophid.utoronto.ca/mirDIP/). Pearson correlation analysis was conducted to screen genes co-expressed with LINC00941 or BASP1-AS1 in the TCGA dataset with r>0.4 and p<0.001. These genes were intersected with predicted mRNAs to define potential mRNA targets of LINC00941 or BASP1-AS1. Finally, the lncRNA-miRNA-mRNA network was constructed using Cytoscape 3.8.0.

Given the fact that there were no widely used LGG cell lines, and glioblastoma and LGG belonged to gliomas, we used two glioblastoma cells (DBTRG-05MG and U251 MG) in this study to explore the effect of HRLs on gliomas. DBTRG was cultured in RPMI-1640 medium with 10% fetal bovine serum (FBS, ExCell Bio, China), and U251 was cultured in DMEM medium with 10% FBS (ExCell Bio, China). Specific siRNAs targeting LINC00941 and BASP1-AS1 were designed and synthesized by GenePharma. Lipofectamine 3000 (Invitrogen, USA) was used for cell transfection.

Total RNA was extracted using TRIzol reagent (Invitrogen, USA). The PrimeScript RT reagent Kit (RR047A, Takara) was used to synthesize cDNA. The TB Green Fast qPCR Mix (RR430S, Takara) was used for qPCR. GAPDH was used as the reference gene. Primers used in this study were as follows: GAPDH: F, CAGGAGGCATTGCTGATGAT; R, GAAGGCTGGGGCTCATTT. LINC00941: F, ACCACTACACTCAGCCAAATAC; R, GGCTATCAACTGTCTCCTTTAGAC. BASP1-AS1: F, AGCACCGGGACACAGAATAG; R, TTTGCGGGAAGGTAAAATTG.

In each well of 96-well plate, 2×10³ cells were inoculated and maintained in culture medium for 0, 24, 48, and 72 h. Then, 10 μl CCK-8 reagent (DOJINDO, Japan) was added into each well, and the optical density of 450 nm was estimated.

In each well of six-well plate, 800 cells were seeded and maintained in culture medium for 14 days. After the fix with 4% paraformaldehyde, cells were stained with 0.1% crystal violet.

A total of 1×10⁵ cells was seeded in 20 mm round coverslip. The EdU Cell Proliferation Assay Kit (Ribobio, China) was used for EdU assay and the processes were conducted according to the manufacturer’s instructions.

Statistical analyses and visualization were mainly performed using R version 3.6.0 and GraphPad Prism version 8.0.1. Student’s t test and one-way ANOVA analysis were used to estimate the differences between two groups and more than two groups. Kaplan-Meier analysis was conducted to compare the survival differences between two groups of patients. Multivariate Cox analysis was used to evaluate the prognostic value of risk score. The correlation of gene expression was determined by Pearson correlation analysis. Time-dependent receiver operating characteristic (ROC) curve analysis was adopted to estimate the predictive value of risk score. Two-sided p ≤ 0.05 was regarded as statistically significant.

To preliminarily explore the role of hypoxia and lncRNAs in LGG, we screened HRLs in the TCGA dataset. A total of 14,488 lncRNAs and 26 hypoxia-related genes were identified and selected for our study ( Figure 1A ). Pearson correlation analysis identified 399 HRLs (|R|>0.5 and p<0.01) in the TCGA dataset ( Table S1 ). After the intersection with lncRNAs in the CGGA dataset, 340 HRLs were identified. With the application of univariate Cox analysis and Lasso analysis, 19 HRLs were selected for further analysis ( Figure 1B ). Among these lncRNAs, four lncRNAs (AL391834.1, LINC00836, BASP1-AS1, and AL023806.1) played protective roles in LGG, whereas the other 15 lncRNAs were risk factors (p<0.05) ( Figure 1C ). Kaplan-Meier analysis further validated the prognostic values of these lncRNAs in the two datasets (p<0.05) ( Figure S1 ). The correlation between 19 selected HRLs and 26 hypoxia-related genes were shown by the heatmap ( Figure 1D ). These results indicated that identified lncRNAs were correlated with hypoxia and were prognostic biomarkers of LGG.

Then we constructed the HRL signature to further characterize the role of HRLs in LGG. With the application of Lasso analysis, the coefficient of 19 HRLs was determined, in which the coefficient of four lncRNAs (AL391834.1, LINC00836, BASP1-AS1, and AL023806.1) was negative and that of the other lncRNAs was positive ( Figure 2A ). The risk score of each LGG patients was calculated according to the coefficient and expression of 19 HRLs. Risk score was significantly associated with the expression of hypoxia-related genes ( Figure S2 ). Then, LGG patients were divided into high-risk and low-risk groups based on the medium value of risk score ( Figures 2B, C ). Patients whose survival time was relatively short and status was censored were enriched in high-risk group ( Figures 2B, C ). The area under curve (AUC) of risk score in predicting 1-, 3-, and 5-year survival of LGG patients was 0.862, 0.874, and 0.805, respectively, in the TCGA dataset, whereas those were 0.835, 0.860, and 0.845, respectively, in the CGGA dataset ( Figure 2D ). Four protective lncRNAs (AL391834.1, LINC00836, BASP1-AS1, and AL023806.1) were highly expressed in low-risk group, whereas other lncRNAs were highly expressed in high-risk group ( Figures 2E, F ). Regarding several well-known biomarkers of glioma, risk score was significantly elevated in grade III glioma compared with grade II one (p<0.05) ( Figure 2G ). In the meantime, risk score was significantly lower in IDH mutant and 1p19q co-deleted gliomas (p<0.05) ( Figures 2H, I ). However, as for the methylation of MGMT promoter, risk score did not share the consistent trend in the TCGA and CCGA dataset, where risk score was significantly elevated in MGMT unmethylated glioma in the TCGA dataset (p<0.05) but not in CCGA dataset (p>0.05) ( Figure 2J ). These findings suggested that the constructed risk signature based on the expression of HRLs was associated with clinical features of LGG and could predict the survival time of LGG patients.

To further verify the prognostic value of risk signature, we constructed a nomogram model, whose C-indexes were 0.857 and 0.772 in the TCGA and CGGA datasets, respectively ( Figure 3A ). Moreover, a calibration plot for probability of survival exhibited satisfactory concordance with the prediction of 3-year and 5-year OS in the TCGA dataset ( Figure 3B ). Then we conducted a subgroup analysis to verify the prognostic value of HRL signature in different subgroups of LGG patients. Results showed that LGG patients with high risk score had poor prognosis (p<0.05) ( Figure 3C ). As for grade II and III gliomas, high risk score indicated worse prognosis (p<0.05) ( Figures 3D, E ). Similarly, in IDH mutant or wildtype LGG patients, those in low-risk group had longer survival time (p<0.05) ( Figures 3F, G ). Moreover, when patients were divided into young (age ≤45 years old) and old (age >45 years old) groups, the prognostic value of risk score was consistent (p<0.05) ( Figures 3H, I ). Furthermore, multivariate Cox analysis revealed that risk score and grade were independent risk factors for LGG patients in the TCGA and CGGA datasets (p<0.05) ( Table 2 ). These results indicated that risk score was a potent marker to predict the prognosis of LGG patients.

Since immunotherapy was a promising therapeutic approach in cancer treatment, we explored the association between risk stratification and the efficacy of ICIs. The expression of several immune checkpoints including PD-1, PD-L1, CTLA-4, TIM-3, B7-H3, IDO1, and LAG3 was significantly elevated in high-risk group compared with the low-risk group in the TCGA and CGGA datasets (p<0.05) ( Figures 4A, B ). Patients with low risk score and low PD-1 expression had significantly better prognosis than those with high risk score and low PD-1 expression (p<0.05) ( Figure 4C ), and patients with low risk score and high PD-1 had prolonged survival than those with high risk score and high PD-1 (p<0.05) ( Figure 4C ). Similarly, the stratification based on HRL signature and immune checkpoints was associated with significant survival difference in LGG patients. Patients with low risk score tended to have better prognosis no matter when the immune checkpoints (PD-L1 and CTLA-4) were highly or lowly expressed (p<0.05) ( Figures 4D, E ). With the application of TIDE algorithm, we found that TIDE score was significantly elevated in high-risk group, which indicated that patients in high-risk group had worse responses to immunotherapy compared with low-risk group (p<0.05) ( Figure 4F ). Thereafter, we extracted the data from IMvigor210 dataset, a cohort of atezolizumab for the treatment of urothelial carcinoma, to investigate the correlation between HRL signature and immunotherapeutic efficacy. Due to the limited number of lncRNA in IMvigor210 dataset, only 11 out of 19 HRLs were identified, and risk score was calculated based on the coefficient and expression of 11 HRLs. Results showed that HRL signature was significantly correlated with the survival of urothelial carcinoma patients receiving atezolizumab treatment (p<0.05) ( Figure 4G ). However, in the IMvigor210 cohort, patients with high risk score had a better prognosis, which might be due to the different roles of HRLs in different types of cancer. Nevertheless, these findings suggested that HRL signature might be a potential biomarker to predict the therapeutic response of immune checkpoint inhibitors.

Further we explored the potential pathways that were associated with the prognosis of patients in high-risk and low-risk groups. GSVA analysis revealed that hypoxia inducible factor 1 (HIF-1) signaling pathway and immune-related pathways including lymphocyte activation, interleukin-mediated signaling pathway, and antigen processing and presentation were highly enriched in high-risk group ( Figures 5A, B ). Differentially expressed genes between high-risk and low-risk groups were screened in the TCGA and CGGA datasets ( Figure S3 ). GSEA analysis showed that genes highly expressed in the high-risk group were enriched in hypoxia-related pathway ( Figure 5C ). Meanwhile, these genes were involved in antigen processing and presentation as well as interleukin secretion in GO pathways ( Figure 5D ); in KEGG terms, these genes were associated with antigen processing and presentation, T cell differentiation, and B cell receptor signaling pathway ( Figure 5E ). Therefore, risk stratification based on risk scores was associated with hypoxia-related pathways, and the high-risk group exhibited highly activated immune characteristics.

Since tumor immune microenvironment was implicated to be associated with the prognosis of LGG patients (25, 26), we also explored the correlation between risk signature and immune microenvironment in LGG. The risk score was significantly positively associated with the stromal score, immune score, and ESTIMATE score (p<0.05) ( Figures 6A–C ). Then tumor purity was calculated according to the algorithm based on ESTIMATE score (21), and high risk score notably indicated low tumor purity (p<0.05) ( Figure 6D ). The infiltration of immune cells was estimated by the conduct of ssGSEA and TIMER algorithms, which contained 28 and six immune cells, respectively. Immune cells such as macrophages, activated T cells, activated B cells, and activated dendritic cells were enriched in high-risk group of samples ( Figures 6E, F ). Besides, risk score was significantly correlated with the abundance of dendritic cell, macrophage, and CD4+ T cells (p<0.05) ( Figures 6G, H ). In addition, CIBERSORT algorithm revealed that the abundance of B cells, macrophages (M1 and M2 subtypes), and naïve CD4+ T cells were highly infiltrated in high-risk group (p<0.05) ( Figures 6I, J ). These results suggested that risk signature was associated with immune infiltration, and high risk score implied low tumor purity, which might account for its risk role in the prognosis of LGG patients.

To investigate the role of HRLs in LGG, we selected two HRLs (LINC00941 and BASP1-AS1), which had the largest absolute coefficient in risk signature, for further analysis. In the TCGA and CGGA datasets, low expression of LINC00941 and high expression of BASP1-AS1 indicated favorable prognosis in LGG patients (p<0.05) ( Figures 7A, B ). Moreover, the expression of LINC00941 was significantly elevated in grade III glioma than grade II one, whereas that of BASP1-AS1 exhibited opposite expression pattern (p<0.05) ( Figures 7C, D ). Similarly, LINC00941 was highly expressed in IDH wildtype glioma, and BASP1-AS1 was highly expressed in IDH mutant glioma (p<0.05) ( Figures 7E, F ). Therefore, LINC00941 and BASP1-AS1 were significantly associated with the prognosis and clinical features of LGG patients. Besides, we found that LINC00941 was significantly negatively correlated with tumor purity (r=−0.19, p<0.05) whereas BASP1-AS1 was positively associated with tumor purity (r=0.51, p<0.05) ( Figure 7G ). LINC00941 was significantly positively associated with the infiltration of B cell, CD8+ T cell, neutrophil, macrophage, and dendritic cell (|r|>0.1, p<0.05), whereas BASP1-AS1 was negatively associated with the infiltration of B cell, CD4+ T cell, neutrophil, macrophage, and dendritic cell (|r|>0.25, p<0.05) ( Figures 7H–M ). Therefore, LINC00941 and BASP1-AS1 were significantly associated with the prognosis and immune infiltration in LGG.

Then we performed in vitro experiments to verify the pathogenic role of LINC00941 and BASP1-AS1 in glioma cells. Three siRNAs were transfected in U251 and DBTRG cells to inhibit the expression of LINC00941 and BASP1-AS1, in which si-LINC00941#1 and si-LINC00941#2 as well as si-BASP1-AS1#1 and si-BASP1-AS1#3 were selected with the relatively high efficiency (p<0.05) ( Figures 8A, B ). The inhibition of LINC00941 significantly reduced the proliferation of glioma cells, whereas the inhibition of BASP1-AS1 significantly promoted their proliferation rates (p<0.05) ( Figures 8C, D ). Colony formation assay indicated that the knockdown of LINC00941 markedly decreased the colony number, whereas the knockdown of BASP1-AS1 exerted reversed effects (p<0.05) ( Figures 8E, F ). Moreover, EdU assay revealed that the proliferation of glioma cells was suppressed by the inhibition of LINC00941 and promoted by the inhibition of BASP1-AS1 ( Figures 8G, H ). Therefore, two prognostic HRLs, LINC00941 and BASP1-AS1, were associated with the proliferation of glioma cells and were potential therapeutic targets for glioma.

The discovery of ceRNA provided a novel insight into the pathogenic role of lncRNA in cancers. Therefore, we applied bioinformatic analyses to explore the potential miRNA and mRNA targets of LINC00941 and BASP1-AS1. Since lncRNA tended to elevate the expression of mRNA by acting as ceRNA, we screened co-expressed mRNAs of LINC00941 and BASP1-AS1 in the TCGA dataset. After the intersect with predicted miRNA and mRNA targets, the ceRNA network of LINC00941 and BASP1-AS1 was constructed ( Figures 9A, B ). Enrichment analysis revealed that LINC00941-targeting mRNAs were associated with cellular senescence and ion channel activity, whereas BASP1-AS1-targeting mRNAs were enriched in membrane potential, synaptic activity, and ion channel activity ( Figures 9C, D ). Differentially expressed genes were screened between high and low expression of LINC00941 or BASP1-AS1 groups ( Figures S3C, D ). GSEA revealed that epithelial mesenchymal transition (EMT), K-ras signaling, reactive oxygen species pathway, and TNF-α signaling pathway were highly enriched in high-LINC00941 group, whereas myc and Wnt/β-catenin were enriched in low-LINC00941 group ( Figure 9E ). Meanwhile, Hedgehog and K-ras signaling pathways were enriched in high-BASP1-AS1 group, and EMT, IL-6/JAK/STAT3, Interferon-α, and TNF-α signaling pathways were enriched in low-BASP1-AS1 group ( Figure 9F ). Therefore, these results indicated that LINC00941 and BASP1-AS1 might affect the proliferation of glioma cells by regulating ion channel activity via modulating EMT and TNF-α signaling pathway.