An Immune Gene-Related Five-lncRNA Signature for to Predict Glioma Prognosis

Background The tumor immune microenvironment is closely related to the malignant progression and treatment resistance of glioma. Long non-coding RNA (lncRNA) plays a regulatory role in this process. We investigated the pathological mechanisms within the glioma microenvironment and potential immunotherapy resistance related to lncRNAs. Method We downloaded datasets derived from glioma patients and analyzed them by hierarchical clustering. Next, we analyzed the immune microenvironment of glioma, related gene expression, and patient survival. Coexpressed lncRNAs were analyzed to generate a model of lncRNAs and immune-related genes. We analyzed the model using survival and Cox regression. Then, univariate, multivariate, receiver operating characteristic (ROC), and principle component analysis (PCA) methods were used to verify the accuracy of the model. Finally, GSEA was used to evaluate which functions and pathways were associated with the differential genes. Results Normal brain tissue maintains a low-medium immune state, and gliomas are clearly divided into three groups (low to high immunity). The stromal, immune, and estimate scores increased along with immunity, while tumor purity decreased. Further, human leukocyte antigen (HLA), programmed cell death-1 (PDL1), T cell immunoglobulin and mucin domain 3 (TIM-3), B7-H3, and cytotoxic T lymphocyte-associated antigen-4 (CTLA4) expression increases concomitantly with immune state, and the patient prognosis worsens. Five immune gene-related lncRNAs (AP001007.1, LBX-AS1, MIR155HG, MAPT-AS1, and LINC00515) were screened to construct risk models. We found that risk scores are related to patient prognosis and clinical characteristics, and are positively correlated with PDL1, TIM-3, and B7-H3 expression. These lncRNAs may regulate the tumor immune microenvironment through cytokine–cytokine receptor interactions, complement, and coagulation cascades, and may promote CD8 + T cell, regulatory T cell, M1 macrophage, and infiltrating neutrophils activity in the high-immunity group. In vitro, the abnormal expression of immune-related lncRNAs and the relationship between risk scores and immune-related indicators (PDL1, CTLA4, CD3, CD8, iNOS) were verified by q-PCR and immunohistochemistry (IHC). Conclusion For the first time, we constructed immune gene-related lncRNA risk models. The risk score may be a new biomarker for tumor immune subtypes and provide molecular targets for glioma immunotherapy.


INTRODUCTION
Glioma is a primary malignant tumor derived from glial cells in the central nervous system. Its annual incidence rate is 7.08 per 100,000 people, and accounts for about 75% of whole brain and other central nervous system malignancies (Lapointe et al., 2018;Ostrom et al., 2019). Clinically, gliomas are often divided into low-grade gliomas (LGGs) and glioblastomas, which have different treatment methods and prognoses. For example, LGGs are slow growing and are mainly treated by total surgical resection. The patient prognosis is relatively good (Sturm et al., 2017). However, the median survival period is less than 2 years with malignant glioblastoma progression, even with standard treatment (surgical resection, adjuvant radiotherapy, and chemotherapy) (Tan et al., 2020). In 2016, the WHO classified gliomas into five categories based on their morphology and molecular characteristics (Louis et al., 2016). Recently, immunotherapy has been used in clinical applications. However, the overall prognosis of glioblastoma patients varies greatly. This may be due to the formation of unique tumor microenvironments during long-term tumor formation and limited molecular markers that distinguish tumor subtypes (Hanahan and Weinberg, 2011). Therefore, it is important to understand the glioma immune microenvironment and screen new molecular markers, which will guide future glioma treatment.
The extracellular matrix, soluble molecules, and tumor stromal cells are the basic components of the tumor microenvironment (Cui et al., 2017). Immune cells and stromal cells are the most common non-tumor cells. Macrophages are the most abundant immune cells in brain tumors (Quail and Joyce, 2017). Glioma often recruits T cells, bone marrow-derived suppressor cells, and macrophages through several pathways to promote immune cell accumulation and transformation into different cell types (Gieryng et al., 2017). Microglia and macrophages are often activated to control anti-tumor immune responses, promote tumor cell proliferation and invasion, and achieve immune escape (Hambardzumyan et al., 2016). Human leukocyte antigen (HLA) (Machulla et al., 2001), programmed cell death-1 (PDL1) (Jackson et al., 2019), cytotoxic T lymphocyte-associated antigen-4 (CTLA4) (Nduom et al., 2015), T cell immunoglobulin and mucin domain 3 (TIM-3) (Das et al., 2017), and other immune-related genes participate in the immune escape process. Therefore, treatments targeting immune checkpoints, microglia, and macrophages are used in the clinic (Poon et al., 2017). However, some patients are in a state of immune tolerance. To improve the quality of medical care and increase the understanding of the immune microenvironment, tumor immune gene analysis is common. Considering tumor-associated immune genes, investigating immune gene sets with guided evolutionary simulated annealing (GESA) can more comprehensively reflect the glioma immune microenvironment in vivo to better establish a prognostic model, find effective molecular markers, and perform effective targeted treatment (Molinaro et al., 2019).
With the development of high-throughput technology and the establishment of public databases, the molecular understanding of tumors has rapidly developed (Serratì et al., 2016), leading to improved understanding of tumor pathogenesis and improved biomarker screening. Importantly, some long non-coding RNA (lncRNA) has been identified as potential glioma biomarkers (Peng et al., 2018). Previously, lncRNAs were hypothesized to have no coding function and were regarded as transcriptional noise. However, lncRNAs play an important regulatory role in gene transcription and posttranscriptional modification. Indeed, lncRNA can regulate inflammation and participate in immune gene expression, thus affecting the tumor immune microenvironment (Chen, 2016;Mathy and Chen, 2017). For example, lincRNA-Cox2 regulates chromatin complex remodeling and participates in inflammatory gene expression (Hu et al., 2016). lncRNA nuclear-enriched abundant transcript 1 (NEAT1) participates in the regulation of interleukin (IL)-8 transcription, thus affecting cytokine response, and induces immune gene expression (Hirose et al., 2014). High HOTAIR lncRNA expression promotes the secretion of monocyte chemoattractant protein-1 (MCP-1/CCL2) by tumor cells and promotes the proliferation of tumor-associated macrophages (TAM) and myeloid-derived suppressor cells (MDSC) in the immune microenvironment (Botti et al., 2019). The complex relationship between lncRNAs and the tumor immune microenvironment has been gradually revealed, and the mechanism of immune-related lncRNA in a variety of tumors has been reported (Hu et al., 2019). However, the relationship between lncRNAs and the glioma immune environment remains unclear.
We analyzed glioma samples downloaded from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA), to examine the glioma immune microenvironment using the single-sample GSEA method. Then, we screened lncRNAs related to the analyzed immune gene set. Using survival curve and Cox regression analysis, a five-lncRNA prognosis model related to the immune gene set was constructed, and the relationship between the risk score and the glioma patient prognosis was explored. Our results provide new ideas for the clinical immunotherapy of glioma.

Patient and Glioma Samples
This study was approved by the patients and the Ethics Committee of the First Affiliated Hospital of Harbin Medical University. All glioma tissue samples were obtained from the surgical resection tissue of glioma patients (n = 18); non-tumor brain tissue was used as the negative control group (n = 5). Tissue samples are stored separately in liquid nitrogen and paraffin embedded.

Immune Grouping and Correlation Analysis
In the single-sample GSEA method, each sample was scored according to 29 immune gene sets and divided into three groups by hierarchical clustering (Molinaro et al., 2019). We used Estimate package to calculate the tumor microenvironment indicators for each sample and analyze the tumor purity (Yoshihara et al., 2013). Then, we used the R-x64-4.0.2 language package to analyze the three immune-related gene and patient prognosis groups. Finally, we analyzed immune cell infiltration in each tumor sample using the CIBERSORT method (Newman et al., 2015) (p < 0.05).

Risk Model
Nine lncRNAs were screened based on the correlation between identified lncRNAs and the immune gene sets (R 2 ≥ 0.62) in CGGA. An additional five prognosis-related lncRNAs were identified using univariate and multivariate survival analyses by Cox regression model (Malinchoc et al., 2000). We divided the samples from the CGGA database into high-and low-risk groups according to the median risk score (Risk score = correlation_lncRNA1 × expression_lncRNA1 + correla-tion_lncRNA2 × expression_lncRNA2 + correlation_lncRNAn × expression_lncRNAn) (Chen et al., 2007;Zhang et al., 2020b). Survival curve and Cox regression analysis were used to construct the immune gene set-related lncRNA risk model.

Risk Model Assessment
We used cor.test function to detect the relationships between lncRNAs (Zhang et al., 2020a). Then, we evaluated the accuracy of the risk model using univariate, multivariate, and receiver operating characteristic (ROC) curves. ggpubr package was used to show the relationship between lncRNAs, clinical symptoms, and immune status. Then, we use principal component analysis for model clustering through scatterplot3d package (Ma and Dai, 2011).

GSEA for Enrichment Analysis
We used clusterProfiler, colorspace, and enrichplot package to perform GO and KEGG analysis based on the sequence of genes which was sorted each gene in descending order of log2FoldChange [log2 (Mean of high immune group genes/Mean of low immune group genes)], and drew a bubble chart (p < 0.05) through ggplot2 package (Cheng et al., 2019).

Quantitative RT-PCR (qRT-PCR)
Total RNA was prepared using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's instructions. The concentration of the total RNA was detected by NanoDrop 2000 (Thermo Scientific TM ). Total RNA (1000 ng) was reverse transcribed into cDNA using qPCR RT Kit (TOYOBO, Japan). Relative expression of target gene to the housekeeping gene GAPDH was determined by qRT-PCR using FastStart Universal 96 SYBR Green Master (ROX) (Roche, Germany). All primer sequence used in this study is listed in Supplementary Table 1. Analysis between the two groups was performed by an unpaired t-test; P < 0.05 was considered statistically significant.

Immunohistochemistry (IHC)
The tissue sample immersed in formalin is wrapped in paraffin and sliced into 5 µm thick sections. Then sample sections were incubated for PDL1, CTLA4, CD3, CD8, and INOS primary antibodies at 4 • C overnight and secondary antibodies at 37 • C for 30 min. Next, samples were visualized by using the diaminobenzidine (DAB) substrate kit for 10 min. After intensive washing, samples were counterstained with hematoxylin, then dehydrated and coverslipped according to manufacturer's protocol. The results of immunohistochemistry (IHC) were taken with Leica microscope.

RESULTS
The Tumor Immune Microenvironment Reflects Tumor Purity Normal brain tissue maintains a low-medium immune state, while gliomas are clearly divided into low-immunity groups (immunity_L), medium-immunity groups (immunity_M), and high-immunity groups (immunity_H) (Supplementary Figure 1A and Figures 1A,B). From immunity_L to immunity_H, the stromal score, immune score, and estimate score (stromal score combined with immune score) increase, and the tumor purity decreases. We further quantified different immunity groups scores and drew violin plots. The changes of immune stromal cells in the tumor microenvironment and the decrease in tumor purity are consistent with Figures 1C,E,  Supplementary Figure 1B, Figure 1G (TCGA, p < 0.001), Figures 1D,F, Supplementary Figure 1C, and Figure 1H (CGGA, p < 0.001). In order to better understand the tumor microenvironment and find potential therapeutic targets, whether there are differences in immune-related genes is worthy of our further study.

Immune Gene Expression in the Three Groups
We generated boxplots to evaluate the expression of immunerelated genes during the immune response. As shown in Figures 2A,B, HLA-related gene expression gradually increased from the immunity_L to immunity_H groups (p < 0.001). We also found that PDL1 (

Patient Prognosis in the Different Immune Groups
To analyze the effect of gene expression on patient prognosis in the different immune groups, we drew survival curves for the TCGA (669 samples: 510 LGG and 159 GBM samples) and the CGGA (971 samples: 596 LGG and 375 GBM samples) (Figure 4). Among the glioma patients, patients in the immune_L group had the best prognosis, followed by the immune_M group, and the immune_H group had the worst prognosis ( Figure 4B, p < 0.001). Among LGGs, prognosis in the different immune groups was similar (Figures 4C,D, p < 0.001). In GBM, the prognosis in the immune_H group tended to be worse than in the immune_L group, but the difference was not statistically significant (p > 0.05).

The Correlation Between lncRNA and Immunity
Using correlation analysis, we found that the lncRNAs in the risk model are associated (Supplementary Figure 3A). The risk score is closely related to the lncRNAs, PDL1, TIM-3, and B7-H3 (Supplementary Figures 3B-I). In addition, we found that AP001007.1, LBX2-AS1, and MIR155HG had the highest expression, while MAPT-AS1 and LINC00515 expression was the lowest in the immune_H group. In contrast, the expression of AP001007.1, LBX-AS1, and MIR155HG were relatively low, while the expressions of MAPT-AS1 and LINC00515 were relatively high in the immune_L group ( Figure 6H). We next analyzed the immune-infiltrating cells in each group. In the immune-H group, we found that naive B cells, plasma cells, CD8 + T cells, regulatory T cells (Tregs), M1 macrophages, M2 macrophages, resting mast cells resting, and infiltrating neutrophils increased. CD4 + naive T cell, inactivated

GO Enrichment and KEGG Pathway Analysis
We used GSEA to analyze enriched differential genes in the immune_H and immune_L groups. We observed that the differential genes were enriched in immunoglobulin complex, circulating, immunoglobulin receptor binding, and MHC protein complex (Figures 7C-F). Further KEGG function analysis (Figure 8A and Supplementary Figure 3J) showed allograft rejection, asthma, intestinal immune network for IgA production, and cytokine-cytokine receptor interaction may be activated cell signaling pathways. The intersection of the two data sets revealed 81 cell signal pathways involved in glioma ( Figure 8B). The five lncRNA we identified may regulate the immune microenvironment through cytokine-cytokine receptor interaction, antigen processing and presentation, complement and coagulation cascades, and intestinal immune network for IgA production (Figures 8C,D).

In vitro Validation of the Risk Model
Through qRT-PCR, we confirmed that AP001007.1, MIR155HG, and LBX2-AS1 are highly expressed, while LINC00515 and MAPT-AS1 are low expressed in gliomas compared to the control group ( Figure 9A). Then, we further found that the risk score was positively correlated with the expression of PDL1, CTLA4, CD3, CD8, and INOS ( Figure 9B, cor > 0.5). Finally, the immunohistochemical results also confirmed that the expression of PDL1, CTLA4, CD3, CD8, and INOS in the high risk group was significantly higher than low risk group base on protein level ( Figure 9C).

DISCUSSION
Glioma cells form a complex regulatory network via the extracellular matrix, stromal cells, and infiltrating immune cells (Hanahan and Coussens, 2012). Some cells secrete factors and lncRNAs to promote inflammation and angiogenesis in tumors, thereby promoting malignant tumor progression and immune escape (Pitt et al., 2016;Dagogo-Jack and Shaw, 2018). It is critical to understand the tumor immune microenvironment and screen new markers to enable targeted glioma therapy for glioma. Abundant research on immune cells has been performed (Charoentong et al., 2017;Jia et al., 2018). However, the cell types, functions, and pathways associated with glioma remain unclear. Therefore, we analyzed 1715 glioma and 1152 normal brain tissue samples using the single-sample GSEA method. We found that the immune environment in gliomas was very different from the immune environment in normal brain tissues. In the immune_H group, the tumor immune cell, and stromal cell content increases, the tumor purity decreases, and tumor heterogeneity becomes greater. These conclusions are in line with previous findings (Hanahan and Coussens, 2012;Pitt et al., 2016;Dagogo-Jack and Shaw, 2018) indicating that this method can accurately reflect the basic conditions of the tumor microenvironment.
Human leukocyte antigens and immune checkpoints are an indispensable regulator of the immune microenvironment (Topalian et al., 2016;Pereira et al., 2019). In the immune_H group, PDL1, CTLA4, TIM-3, and CD96 expression were increased. Immune checkpoints act to negatively regulate immune regulation. Normal immune surveillance and cell killing ability are weakened in many tumors. Further, tumors often have immune escape or immunotherapy resistance mechanisms, leading to ineffective clinical treatment (Field et al., 2017;Qian et al., 2018). We also observed that among all the glioma samples, the immune_H group had the worst prognosis, followed by the immune_M group, and the immune_L group had the best prognosis. This conclusion also supports previous results. Moreover, GSEA has been used in many studies and has a certain degree of credibility (Ma et al., 2020;Wang et al., 2020). Therefore, the single-sample GSEA method based on expressed immune genes can distinguish the biological characteristics of the immune microenvironment between different gliomas, which provides the possibility for screening the immune microenvironment-related biomarkers.
Long non-coding RNAs can regulate the tumor immune microenvironment and can be used as biomarkers. For example, NF-kappaB interacting lncRNA (NKILA) can promote the immune escape of tumor cells by regulating T cell activity (Huang et al., 2018). Further, SATB2-AS (the antisense transcript   of SATB2 protein) can directly combine with WDR5 (WD repeat containing protein 5) and GADD45A (growth arrest and DNA damage protein 45A) to regulate SATB2 expression, thereby inhibiting tumor cell metastasis and regulating the tumor immune microenvironment . Immune-related lncRNAs have carcinogenic effects in several tumors, and can be used as biomarkers . Thus, it is very important to determine whether lncRNAs related to the immune gene set have clinical diagnostic and prognostic value. We screened five lncRNAs using the Cox regression method and constructed a prognostic model. We found that the risk score is related to prognosis and is an independent factor that can be used for clinical diagnosis. We further observed that five lncRNAs interact and are closely related to the clinical symptoms of glioma patients (WHO grade, IDH1 status, 1q19q status, and MGMT). Principle component analysis (PCA) analysis showed that subgroups within the high-and low-risk groups can be well distinguished using our method. These conclusions show that the risk scores of the five lncRNAs related to the immune gene set can predict patient prognosis and clinical characteristics, and can be used as a new biomarker to inform clinical diagnosis and treatment. We used boxplots to visualize lncRNA expression in each immune group. In the immune_H group, we found that AP001007.1, LBX2-AS1, and MIR155HG were highly expressed, while LINC00515 and MAPT-AS1 expression was low. The immune_L group showed the opposite trend. Survival analysis showed that AP001007.1, LBX2-AS1, and MIR155HG were risk factors, and their high expression predicted poor patient outcomes. LINC00515 and MAPT-AS expression were protective indicators, and low expression predicted poor patient prognosis. LBX2-AS1 produces malignant behavior in gliomas by conferring resistance to cell apoptosis (Chen et al., 2020). MIR155HG promotes immune cell infiltration and immune resistance (Peng et al., 2019). In contrast, MAPT-AS1 expression indicates a good prognosis for cancer patients . Therefore, the immune gene set-related model we constructed has considerable credibility.
We found that the five lncRNAs we analyzed may promote immune cell infiltration through cytokine-cytokine receptor interaction, antigen processing and presentation, complement and coagulation cascades, and may contribute to immune resistance and tolerance, ultimately leading to poor patient prognosis. However, this study also has some limitations. First, because there was no information regarding MGMT, 1p19q, and other related molecules in the TCGA dataset, only a single CGGA cohort was used for statistical analysis. The CGGA sample information is mainly clinical sample information from Chinese patients, which may only indicate region-specific effects. Second, basic experiments have also verified the important role of some immune gene-related lncRNAs in regulating glioma development. However, the mechanism underlying our prognostic model remains unclear and requires additional experiments to verify our in silico results. Third, we confirmed that the five lncRNA have potential clinical value to identify risk factors, but more factors should be considered, especially considering multimodal glioma development.

CONCLUSION
Immune-related genes can reflect the characteristics of the immune microenvironment. To reveal the mechanism of partial resistance or treatment resistance within a new risk model, five immune-related lncRNAs were analyzed and shown to have good stability and feasibility (AP001007.1, LBX-AS1, MIR155HG, MAPT-AS1, and LINC00515). Thus, our study reveals biomarkers that distinguish specific glioma groups and can be used in the clinical diagnosis and treatment of glioma.

DATA AVAILABILITY STATEMENT
The public databases mined in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
This study was approved by the Ethics Committee of the First Affiliated Hospital of Harbin Medical University.

AUTHOR CONTRIBUTIONS
XZW, MG, JYY, and QY collated and analyzed the data. CZ, QYJ, STW, DYH, XH, and LGW completed the writing and repair of the manuscript. JZ, HZ, JNW, and SGZ designed and guided the subject. All authors reviewed and approved the final manuscript.