Comprehensive Analysis of Genome-instability Related lncRNAs and Tumor-immune Microenvironment to Discriminate the Prognosis, Immunotherapy and PARP Inhibitor Responses in Patients with Low-grade Glioma

Background: Previous study revealed that Genome-instability was correlated with tumor-immune microenvironment in cancer. We try to discriminate the prognosis, immunotherapy and poly (ADP)-ribose polymerase (PARP) inhibitor responses through comprehensive analysis of genome-instability related lncRNAs and the tumor-immune microenvironment in patients with low-grade glioma (LGG). Methods and Results: RNAseq data, genome variation proling data and copy number variation (CNV) data were used to evaluate the genomic instability of LGG patients. Genomic unstable-like (GU-like) and genomic stable-like (GS-like) clusters were identied by hierarchical clustering analysis of 102 genome-instability related lncRNAs (GILncRNAs). GS-like cluster had a tendency to receive better clinical outcome. Patients in GU-like cluster were more likely to respond to immunotherapy, especially anti-PD-1/PD-L1 treatment. PARP inhibitors including Rucaparib and Olaparib will get better therapeutic effects for patients in GU-like cluster. Lasso and Cox regression analysis were utilized to construct the risk model based on GILncRNAs. As for the risk model constructed by 9 GILncRNAs, the overall survival, clinical outcome, immunotherapeutic response, and PARP inhibitor sensitivity were signicantly different between patients of high and low-risk groups. Conclusions: The genome-instability related lncRNAs signature involved in our risk model had great advantages in predicting prognosis, immunotherapy and PARP inhibitor response. Univariate Cox regression analysis was rstly employed to evaluate the correlation between expression level of genome-instability related lncRNAs (GILncRNA) and overall survival of patients with LGG using “survival” R package, in which COX p value less than 1.0*10^ -9 was considered as lncRNAs with prognostic value. The results were visualized by “forestplot” R package. The LGG samples were randomly divided into train and test group at a ratio of 1:1. The least absolute shrinkage and selection operator (LASSO) regression algorithm with penalty parameter tuning conducted by 1000-fold-crosss validation was used to further select the valuable GILncRNAs from the previous obtained prognostic lncRNAs. The Receiver Operating Characteristic (ROC) curve analysis was conducted and the Areas Under Curve (AUC) values were calculated to evaluate the accuracy of the risk model. R packages including "survival","caret","glmnet","survminer" and "timeROC" were involved above. The risk score for each sample was calculated by the following formula: This study provided a comprehensive analysis of genomic instability and tumor-immune microenvironment in low-grade glioma through exploring the expression patterns of genome-instability related long non-coding RNAs (GILncRNAs), looking forward to nding a novel method for classication to distinguish subtypes with different prognosis, immunotherapy and PARP inhibitor response. Our results suggested that the GILncRNAs were signicantly related to prognosis, immunotherapy and PARP inhibitor response of LGG patients than other established lncRNA signatures and had great advantages for prediction. Firstly, LGG samples of genomic unstable (GU) group and genomic stable (GS) group were extracted from TCGA dataset, followed which the genomic instability of the two groups were evaluated from the aspect of tumor mutation burdens, mutation frequencies of specic genes and copy number variations. Through the comparison of expression proles of lncRNAs between GS and GU groups, we obtained genome-instability related lncRNAs (GILncRNAs), based on which 529 LGG samples were divided into two clusters, namely GS-like cluster and GU-like cluster. Then, we compared their clinicopathological parameters to clarify the correlation between the clustering and clinicopathological features. Distinct tumor-immune microenvironment (TIME) and immunogenomic patterns were identied between patients of GS-like and GU-like clusters. Patients in GU-like cluster were more likely to respond to immunotherapy, especially anti-PD-1/PD-L1 treatment, as predicted by unsupervised subclass mapping analysis. PARP inhibitors including Rucaparib and Olaparib will get better therapeutic effects for patients in GU-like cluster. In order to facilitate the clinical practice, genome-instability related lncRNAs signature was established as a powerful risk model for predicting the prognosis of LGG patients. According to the calculated risk score, patients were further divided


Introduction
Low-grade gliomas, composed of grade II-III according to the World Health Organization classi cation system, are the most common primary malignant tumors in central nervous system in children, and differ with high-grade glioma (glioblastoma, GBM) in biological and clinical features [1,2]. The therapeutic outcome is not satisfactory for LGG patients received conventional treatment [3,4]. Research focusing on LGG treatment still has a long way to go. Genomic instability is a basic and common feature in the development of malignant tumors, which were characterized as either chromosomal instability (CIN) or microsatellite instability (MSI) [5]. Due to rapid DNA replication and transcription, tumor cells are sensitive to DNA damage, especially damage from double-strand breaks (DSBs) [6]. In response to DNA damage, the DNA damage response (DDR) network was activated to be implicated in the repair of DNA damage, in which erroneous repair is correlated with genomic instability [7,8]. Poly (ADP)-ribose polymerase (PARP) plays a critical role for DNA repair in cancer cells with DDR gene mutations [9] and inhibition of PARP as a synthetic lethal treatment has got promising therapeutic effect in breast cancer [10]. Previous studies have demonstrated an underlying correlation between genomic alterations and tumor-immune microenvironment (TIME) [11,12]. Moreover, TIME has been reported to play a crucial role in the progression and aggressiveness of cancer. The abundance and type of in ltrating immune cell were correlated with cancer therapy and clinical outcome [13,14]. In recent years, an increasing number of studies have been paying attention to immunotherapy such as immune checkpoint inhibitors, antitumor vaccines and cellular immunotherapy for glioma. However, the clinical outcome for patients received immunotherapy can be diverse and only a minority of patients bene tted from it [15,16]. Considering the many factors affecting the response to immunotherapy, a concise method for physicians to predict the immunotherapeutic response for individuals is on need. Long non-coding RNAs (lncRNAs) have been demonstrated to be involved in various biological process as they directly interact with the cell cycle, proliferation pathways and microbiome balance. Furthermore, lncRNAs are increasingly deemed to be important in the regulation of immune cell differentiation and function [17]. A large number of prognostic models based on lncRNAs have been established to predict the prognosis of patients with cancer [18][19][20][21]. Sun et al. identi ed that lncRNAs were related to tumor immune in ltration and built tumor immune in ltration-associated lncRNAs signature for improving prognosis and immunotherapy response of patients with non-small cell lung cancer [22]. However, there is rare study focusing on investigating a risk model for predicting prognosis, immunotherapy and chemotherapy response for patients with LGG from the respect of genome-instability related lncRNAs. In this manuscript, we aimed to construct a precise risk model to predict prognosis, immunotherapy and PARP inhibitor responses through comprehensive analysis of genome-instability related lncRNAs (GILncRNAs) and tumor-immune microenvironment in low-grade glioma.

Screening of lncRNAs Related to Genome Instability
The RNA annotation le of Genome Reference Consortium Human Build 38 (GRCh38) was downloaded from the Ensembl website (http://asia.ensembl.org/) for annotation of RNAseq data. The cumulative counts of somatic mutations in each sample were calculated and ranked in a decreasing manner, followed which the top 25% of the patients and the bottom 25% of the patients were designed as genomic unstable (GU) group and genomic stable (GS) group, separately. The mutation types and frequencies of the top 20 genes with the highest mutation frequency were analyzed and visualized using "maftools" package in R (version 4.1.0). To better exhibit the gain/loss alterations of chromosomes, CNV summary plots were visualized by Circos plots using the "RCircos" R package. LncRNA and mRNA expression data were all log2 transformed for differentially expression analysis using R project. The Page 4/25 differentially expressed lncRNAs between the two groups were analyzed using the Wilcoxon rank-sum test in "limma" package of R software and p values were regulated by false discovery rate (FDR).

Unsupervised Hierarchical Clustering Analysis
Unsupervised hierarchical clustering analysis was conducted through "sparcl", "pheatmap" and "limma" packages of R software using the differentially expressed lncRNAs between GU/GS groups and samples were divided into two gene clusters. The gene cluster with higher mutation counts was identi ed as genomic unstable-like (GU-like) cluster, whereas the other was de ned as genomic stable-like (GS-like) cluster.

Differential Analysis between GU-like and GS-like Clusters
We identi ed the differentially expressed mRNAs between the two gene clusters (|Fold Change| >4.0 and false discovery rate (FDR) adjusted P<0.01]. Top 20 differentially expressed genes (DEGs) were visualized in a heatmap. Subsequently, functional annotation and pathway enrichment analyses, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, were carried out on the DEGs by using "clusterPro ler," "org.Hs.eg.db," "enrichplot," and "ggplot2" packages of R software. FDR <0.001 was considered statistically signi cant in GO analysis and FDR <0.05 was considered signi cantly enriched in KEGG analysis.

Construction of Protein-Protein Interaction (PPI) Network
The PPI of DEGs between GU-like and GS-like clusters was explored by the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database online (https://string-db.org/). The PPI network was conducted by minimum required interaction score ≥ 0:900 with active interaction sources including textmining, experiment, databases, co-expression, neighborhood, gene fusion and co-occurrence. Top 30 key genes were identi ed according to the number of adjacent nodes.

Analysis of Tumor-Immune Microenvironment
The in ltration of 23 immune cells in tumor was estimated through single sample gene set enrichment analysis (ssGSEA) by using "GSEABase" and "GSVA" R package. We obtained the ssGSEA score for each sample according to the given immune gene set le, after which the scores were normalized for further analysis. The abundance of in ltrating immune cells and related immune functions between subgroups were evaluated and visualized in boxplot and heatmap. The correlation between risk score and in ltrating immune cells was analyzed and visualized using "limma" and "corrplot" R packages. Stromal and immune scores were calculated by the "ESTIMATE" R package based on the gene expression pro les of Page 5/25 mapping (https://cloud.genepattern.org/gp/) method was used to predict the response to immunotherapy, such as anti-PD1 and anti-CTLA4 therapy, in the different subtypes.
Construction and Veri cation of Prognostic Signature based on Genome-Instability related LncRNAs Univariate Cox regression analysis was rstly employed to evaluate the correlation between expression level of genome-instability related lncRNAs (GILncRNA) and overall survival of patients with LGG using "survival" R package, in which COX p value less than 1.0*10^-9 was considered as lncRNAs with prognostic value. The results were visualized by "forestplot" R package. The LGG samples were randomly divided into train and test group at a ratio of 1:1. The least absolute shrinkage and selection operator (LASSO) regression algorithm with penalty parameter tuning conducted by 1000-fold-crosss validation was used to further select the valuable GILncRNAs from the previous obtained prognostic lncRNAs. The Receiver Operating Characteristic (ROC) curve analysis was conducted and the Areas Under Curve (AUC) values were calculated to evaluate the accuracy of the risk model. R packages including "survival","caret","glmnet","survminer" and "timeROC" were involved above. The risk score for each sample was calculated by the following formula: where the coefIncRNAi represents the coe cient of ith lncRNA, and the EXP IncRNAi represents the expression level of the ith lncRNA in the risk model. Univariate and multivariate Cox regression analysis were used to examine the independent prognostic value of risk score computed by the model. The distribution of clinical features among high and low risk groups was compared and visualized using "limma" and "ComplexHeatmap" packages in R. In addition, the drug sensitivity was predicted between the subtypes using "pRRophetic" R package, which was determined by the half maximal inhibitory concentration (IC50) of each sample based on the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). Nomogram integrating risk score and clinical factors was utilized to optimize the prognostic model by using "rms" R package. Prediction for 1-, 3-and 5-years survival was identi ed based on multivariate Cox regression results. Calibration curves were used to determine the performance of the nomogram. PCA and t-SNE analysis were employed to examine whether the 9 genome-instability related lncRNAs involved in the risk model can distinguish the samples in our research by using "Rtsne" R package.

Statistical Analyses
The Kruskal-Wallis test was used to compare multi-groups, and the Wilcoxon test was used to compare two groups. Chi-square tests were implemented to evaluate the distribution of categorical variables between subgroups. The student's t-test was employed to compare the continuous data between pairs of subgroups. Correlations between normally distributed variables were analyzed by Pearson's correlation test, while correlations between non-normally distributed variables were assessed with Spearman's correlation test. The Kaplan-Meier plotter was employed to generate survival curves for the subgroups, and the log rank test was used to evaluate the statistically signi cant differences. Two-tailed p < 0.05 was considered as statistically signi cant unless speci cally noted. R version 4.1.0 (Institute for Statistics and Mathematics, Vienna, Austria) was used to conduct the statistical analysis and visualization of the results.

Screening and identi cation of genome-instability related lncRNAs in LGG patients
The clinical and pathological characteristics of the LGG cohort used in our study were displayed in Table 1.Patients were divided into genomic unstable (GU) group and genomic stable (GS) group based on the cumulative number of somatic mutations. In order to evaluate the difference of genomic instability, we investigated the mutation types, mutation frequencies and copy number variations of the two groups.
Top 20 genes with the highest alteration frequencies were analyzed ( Fig. 1a and b). The results of the oncoplots and bar plots revealed that the alteration frequencies of most of the genes were signi cantly higher in patients of GU group than those of GS group including EGFR, FLG, HMCN1, LRP2, NOTCH1, OBSCN, PTEN, RYR2, TTN (The bar plots for FLG, HMCN1, LRP2, NOTCH1, OBSCN, RYR2 and TTN were displayed in Online Resource 5A, P < 0.05). However, the proportion of patients with IDH1 mutations was signi cantly higher in GS group (80%) compared with GU group (54%, p < 0.001) suggesting that IDH1 may play a crucial role in the tumorigenesis and work synergistically with other genes involved in the same pathway [24]. The analysis of CNVs of chromosomes was further performed to explore the distinct genomic variations between the two groups ( Fig. 1c and d). Patients in GU group tended to bear a greater burden of copy number ampli cations (P = 8.9 × e − 13 ) and deletions (P = 3.5×e − 7 ) than those in GS group. All these ndings indicated a signi cant difference in the genomic instability between the two groups. We compared the lncRNA expression pro les between GU and GS groups. 102 genome-instabilityrelated lncRNAs were identi ed with |Fold Change| >2.0 and adjusted p < 0.05. Heatmap of top 20 differentially expressed lncRNAs was shown in Fig. 1e. Table 1 Clinical information for LGG patients in this study. Identi cation of two LGG clusters based on genomeinstability related lncRNAs with distinct survival outcomes and clinical features As shown in the Fig. 2a, unsupervised hierarchical clustering analysis was conducted for LGG patients based on 102 genome-instability related lncRNAs (GILncRNAs). The somatic mutation counts were signi cantly different between the two clusters (Fig. 2b, p < 2.22×e − 16 ). Thus, one cluster with higher mutation counts was de ned as GU-like cluster, the other was de ned as GS-like cluster. The expression level of UBQLN4, which is correlated with genomic instability, was signi cantly higher in patients of GSlike cluster than those of GU-like cluster (Fig. 2b, p = 1.4×e − 13 ). The overall survival (OS) analysis was employed between the two clusters, indicating that patients in GS-like cluster displayed better prognosis compared with GU-like cluster (Fig. 2c, p < 0.001). Heatmap of multi-clinicopathological characteristics was adopted to demonstrate the correlation between the clustering and clinical features (Fig. 2d). As shown in Fig. 2e, patients in GU-like cluster were signi cantly older than that in GS-like cluster (p < 0.001).
Tumors in patients of GU-like cluster were more likely to recur or progress compared with GS-like cluster (p < 0.001). Histologic grade for patients in GU-like cluster was higher than that in GS-like cluster (p < 0.001). Patients in GU-like cluster tended to receive postoperative radiotherapy (p < 0.001). IDH1 for patients in GU-like cluster was more frequently mutated compared with GS-like cluster (p < 0.001). Karnofsky Performance Scores for patients in GU-like cluster were higher than that in GS-like cluster (p = 0.002), indicating self-care abilities for patients in GU-like cluster were more likely to be affected. Clinical outcomes for patients in GU-like cluster seemed to be worse than those in GS-like cluster (p = 0.01), with more patients in tumor-bearing status. Differentially expressed mRNAs (DEGs) were identi ed between the two clusters (|Fold Change| >4.0, adjusted p < 0.01). Heatmap exhibiting top 20 DEGs was shown in Fig. 2f. Protein-Protein interaction network was employed to further evaluate the interplay of the DEGs (Fig. 2g). 30 hub genes with the highest number of adjacent nodes were noted in Fig. 2h.

Enrichment analysis identi ed molecular mechanisms affected by DEGs between GS-like and GU-like clusters
To explore the underlying different molecular mechanisms between the two clusters, GO analysis was employed, regarding molecular function (MF), cellular component (CC), and biological process (BP). As shown in Online Resource 1A-C, the DEGs were signi cantly enriched in molecular functions associated with formation of genomic instability, including DNA-binding transcription activator activity, RNA polymerase II-speci c. KEGG pathway analysis identi ed signaling pathways associated with tumorigenesis and tumor immunology including cytokine-cytokine receptor interaction pathway, ECMreceptor interaction pathway and transcriptional misregulation in cancer (Online Resource 1D-F). These results of functional enrichment analysis suggested that GU-like cluster differs in genome-instability and tumor-immunology associated functions compared with GS-like cluster, implying signi cant differences in tumor-immune environment and genomic instability.
Distinct tumor-immune microenvironment (TIME) and immunogenomic patterns were identi ed between GS-like and GU-like clusters Based on the ESTIMATE algorithm, the compositions of the TIME for each sample were analyzed.
Compared with those of GS-like cluster, the immune scores (p = 2.9xe − 13 ), stromal scores (p < 2.22 x e − 16 ) and ESTIMATE scores (p = 3.3 x e − 15 ) were signi cantly higher in patients of GU-like cluster, indicating higher abundances of immune and stromal cells and lower tumor purity (Fig. 3a). Furthermore, comparisons of the abundance of 23 tumor in ltrating immune cells in different clusters were carried out using ssGSEA. Most of the immune cells were signi cantly more abundant in patients of GU-like cluster (Fig. 3b). The heatmap of the in ltration levels of the immune cells was displayed in the Fig. 3c. As shown in Fig. 3d, the associated immune functions were more active in patients of GU-like cluster.
Patients of GU-like cluster were more sensitive to immunotherapy and PARP inhibitors An increasing number of researches revealed that immunotherapy targeting immune check-points brought hope for clinical treatment of malignant tumors [25]. Identi cation of genome-instability related lncRNA signature for prognosis Univariate cox proportional risk regression was utilized to select the genome-instability related lncRNAs with prognostic values, combining the expression levels of genome-instability related lncRNAs and survival information of LGG patients. 33 genome-instability-related lncRNAs with predictive prognosis roles were determined with p value < 1.0 x 10 − 9 . The results were shown in the forest plot with hazard ratio (HR) of 95% con dence intervals (Fig. 4a). Patients with LGG were randomly divided into train and test group at a ratio of 1:1. Lasso regression analysis was applied to further evaluate and lter these lncRNAs in the train set and 9 lncRNAs were identi ed as the key factors in the risk model with 1000 replicates of cross validation (Fig. 4b and c), including LINC01831, CRNDE, AC025171.5, AL390755.1, AC012073.1, AL355974.2, AC010273.3, AC131097.3, AC092718.3. LncRNSs with associated coe cients were listed in Fig. 4d. The risk score was calculated by the following formula: LGG patients were further divided into high-risk group and low-risk group with the median risk score as the cutoff, in the train and test set, respectively.
Evaluation of genome-instability related lncRNA signature for prognosis Kaplan-Meier analysis showed that overall survival of patients in low-risk group were signi cantly higher compared with those in high-risk group in the train and test set ( Fig. 4e and l, p < 0.001). ROC curve analysis of the risk model based on genome-instability related lncRNAs signature (GILncSig) over time for 1, 2, and 3 years in the train set showed AUC values with 0.892, 0.912, and 0.905, respectively (Fig. 4f). In the test set, AUC values over time for 1, 2, 3 years were 0.819, 0.899 and 0.900, separately (Fig. 4m). In the train set, the AUC value of GILncSig over time for 1 year is 0.870, which is higher than that of other clinicopathological factors including diagnosis age (AUC = 0.704), sex (AUC = 0.509), histologic type (AUC = 0.354) and histologic grade (AUC = 0.669). Similar results were obtained regarding the ROC curve analysis over time of 2, 3 and 5 years (Online Resource 5B). The expression levels of 9 genomeinstability-related lncRNAs in the risk model were signi cantly higher in patients of high-risk group than those of low-risk group ( Fig. 4g and n). The scatter plots of survival time showed that the survival time decreased with the increasing of the risk score and the mortality rate was higher in patients of high-risk group ( Fig. 4h and o). As for the scatter plots of somatic mutations, it seemed that the somatic mutation counts were also higher in patients of high-risk group ( Fig. 4i and p) and it was validated in the boxplot with p value less than 0.01 ( Fig. 4k and r). The expression level of UBQLN4 was signi cantly higher in patients of low-risk group compared with high-risk group (Fig. 4l and s). Univariate and multivariate cox regression were performed to analyze whether the risk score calculated based on the above 9 genomeinstability related lncRNAs were independent prognostic factors of clinical or pathological factors (age, sex, histologic grade, histologic type, and mutation count Validation of genome-instability related lncRNA signature for prognosis GBM dataset from TCGA database was used in this study to verify the performance of the risk model based on 9 genome-instability related lncRNAs. According to the expression levels of the 9 genomeinstability related lncRNAs in GBM dataset, risk score for each sample was calculated, after which patients were divided into high and low-risk group according to the cutoff of median score. Kaplan-Meier curves showed that patients in the high-risk group had signi cantly lower overall survival probability compared with the low-risk group (Online Resource 2G, p < 0.001). As the standard chemotherapy for GBM postoperatively, temozolomide (TMZ) can signi cantly increase the overall survival of GBM patients according to the previous studies. The estimated IC50 values of TMZ were signi cantly lower in patients of low-risk group (p = 0.0071), indicating that GBM patients in high-risk group tended to be more resistant to TMZ therapy than those in low-risk group (Online Resource 2H). There was no statistical signi cance in the estimated IC50 values of 3 PARP inhibitors between patients in low-and high-risk groups (Online Resource 5C).
In addition, the predictive performance of our risk model based on genome-instability related lncRNAs signature (GILncSig) was compared with other glioma-associated prognostic lncRNAs signatures established in previous studies (Online Resource 2I). Although the predictive power of our GILncSig for 1 year rate was not better than Qiu's model, the AUC values of our GILncSig were signi cantly higher than other three risk models over time for 2, 3 and 5 years, indicating that the risk model created by our group had better prognostic performance in predicting overall survival of LGG patients [18][19][20].
Distinct clinicopathological features, tumor-immune microenvironment (TIME), immunotherapy response and sensitivity to PARP inhibitors were identi ed between different risk groups PCA and t-SNE analysis revealed that the 9 genome-instability related lncRNAs involved in our risk model can distinguish the samples with LGG (Online Resource 3A). As shown in the heatmap (Online Resource 3B), the distribution of diagnosis age, disease free time, disease free status, histologic grade, histologic type, IDH1 mutation status, Karnofsky performance score, postoperative radiotherapy status, person neoplasm status, and therapy outcome between the two groups were signi cantly different. The detailed results were demonstrated in boxplots to show the distribution of risk scores of patients with multi-clinicopathological characteristics, indicating that the risk scores of patients with worse clinicopathological features tended to be signi cantly higher (Online Resource 3C).
The distribution of patients in the two gene-clusters corresponding with the risk groups and survival status was shown in alluvial diagram (Online Resource 4A). Comparison between patients of high-and low-risk groups demonstrated that the immune scores (p = 8 x e − 5 ), stromal scores (p = 0.00018) and ESTIMATE scores (p = 8.2 x e − 5 ) were signi cantly higher in patients of high-risk group, indicating higher abundance of immune and stromal cells, while lower tumor purity (Online Resource 4B). Furthermore, most of the immune cells were more abundant in patients of high-risk group (Online Resource 4C and D). The associated immune functions were more active in patients of high-risk group (Online Resource 4E).
The expression levels of all the immune checkpoints mentioned above were elevated in patients of highrisk group compared with low-risk group (Online Resource 4F). In addition, we found a positive correlation between the expression of these immune check-points and risk score (Online Resource 4G). The results of unsupervised subclass mapping analysis demonstrated that patients in high-risk group were more likely to respond to the immunotherapy of anti-PD1 (Fig. 5a, Bonferroni corrected p = 0.007). However, no positive results were found regarding the response to anti-CTLA4 therapy between the two groups. The estimated IC50 values of Rucaparib, Olaparib and Veliparib were signi cantly lower in patients of highrisk group (Fig. 5b-d, p = 1.4×e − 6 , p = 8.3×e − 9 , p = 0.0034, respectively), suggesting that patients in low-risk group tended to be more resistant to PARP inhibitor therapy than those in high-risk group.

Discussion
As a common subtype of glioma, LGG which was composed of astrocytoma, oligoastrocytoma and oligodendroglioma, was not sensitive to current therapy [3]. Accumulating studies have placed special emphasis on the investigation of immunotherapy for glioma, which was revolutionizing cancer care in recent years [26,27]. However, only a minority of patients can get a satisfactory therapeutic effect [28]. It is extremely important to select proper patients for immunotherapy. Taking anti-PD1/PDL1 treatment as an example, there are numerous factors affecting the response to immunotherapy including the expression level of PDL1, tumor mutation burden and abundance of in ltrating immune cells, which makes the clinical decision complex [29,30]. There is an urgent need for an easy method to select proper potential candidates to receive immunotherapy. Recent researches had shown that a large number of long noncoding RNAs (lncNRAs) had important roles on the prediction of prognosis and immune microenvironment of gliomas [18][19][20]22]. So far, there is no ideal classi cation or prognostic model to identify LGG patients with different response to immunotherapy and chemotherapy based on genomeinstability related lncRNAs. Inhibition of PARP is a new strategy of chemotherapy for glioma patients [31]. Inhibition of PARP triggers a synthetically lethal effect in tumors with homozygous deletions, deleterious mutations, or both in DNA damage repair genes including BRCA1/2, ATM, PALB2, and the FANC gene family, which is correlated with genomic instability [32][33][34][35]. PARP inhibitors, such as Rucaparib, Olaparib and Veliparib, can inhibit the activity of PARP so as to prevent DNA replication and transcription as well as generate double-strand DNA breaks, and eventually led to cell death [36]. UBQLN4, a member of the ubiquilin family, which was correlated with genome instability, had been demonstrated to be associated with actionable poly (ADP)-ribose polymerase (PARP) inhibitor sensitivity [37]. Drug sensitivity regarding PARP inhibitors for LGG patients was involved in our study.
This study provided a comprehensive analysis of genomic instability and tumor-immune microenvironment in low-grade glioma through exploring the expression patterns of genome-instability related long non-coding RNAs (GILncRNAs), looking forward to nding a novel method for classi cation to distinguish subtypes with different prognosis, immunotherapy and PARP inhibitor response. Our results suggested that the GILncRNAs were signi cantly related to prognosis, immunotherapy and PARP inhibitor response of LGG patients than other established lncRNA signatures and had great advantages for prediction. Firstly, LGG samples of genomic unstable (GU) group and genomic stable (GS) group were extracted from TCGA dataset, followed which the genomic instability of the two groups were evaluated from the aspect of tumor mutation burdens, mutation frequencies of speci c genes and copy number variations. Through the comparison of expression pro les of lncRNAs between GS and GU groups, we obtained genome-instability related lncRNAs (GILncRNAs), based on which 529 LGG samples were divided into two clusters, namely GS-like cluster and GU-like cluster. Then, we compared their clinicopathological parameters to clarify the correlation between the clustering and clinicopathological features. Distinct tumor-immune microenvironment (TIME) and immunogenomic patterns were identi ed between patients of GS-like and GU-like clusters. Patients in GU-like cluster were more likely to respond to immunotherapy, especially anti-PD-1/PD-L1 treatment, as predicted by unsupervised subclass mapping analysis. PARP inhibitors including Rucaparib and Olaparib will get better therapeutic effects for patients in GU-like cluster. In order to facilitate the clinical practice, genome-instability related lncRNAs signature was established as a powerful risk model for predicting the prognosis of LGG patients. According to the calculated risk score, patients were further divided into high-and low-risk groups. Patients in high-risk group tended to have worse overall survival, while immunotherapy and PARP inhibitors will get better therapeutic effects. In addition, the risk model was also applied to GBM patients to examine the performance in high-grade glioma. Unexpectedly, the overall survival of GBM patients with higher risk scores was signi cantly worse than those with lower risk scores. Moreover, TMZ will get better therapeutic effect for patients with lower risk scores.
In this study, LGG patients were separated into two clusters based on 102 genome-instability related lncRNAs (GILncRNAs). Patients of GS-like cluster had a tendency to receive better prognosis and clinical outcome. On the other hand, the TIME and immunogenomic patterns including stromal score, immune score, ESTIMATE score, abundance of in ltrating immune cells, immune functions and expression of immune checkpoints were signi cantly higher in patients of GU-like cluster, laying foundation for the response to immunotherapy. The expression level of UBQLN4 in patients of GS-like cluster was signi cantly higher than those of GU-like cluster and the somatic mutation counts which can be an indicator of genomic instability were signi cantly lower in patients of GS-like cluster, which were consistent with the results of Shiloh's group. However, the estimated IC50 values of PARP inhibitors including Rucaparib and Olaparib in patients of GS-like cluster were signi cantly higher compared with GU-like cluster, which indicated that patients of GS-like cluster were more resistant to PARP inhibitor treatment. However, Shiloh's group demonstrated that UBQLN4 overexpression in tumors promotes PARP1 inhibitor sensitivity, which seemed differing with our results [37]. Considering they focused on the neuroblastoma cases, the difference of subtype of glioma may account for these results. In addition, indepth investigation on this issue was demanded to clarify the relationship between UBQLN4 and PARP inhibitor sensitivity.
As for the risk model constructed by 9 genome-instability related lncRNAs (GILncRNAs), the overall survival, immunotherapeutic response, PARP inhibitor sensitivity and other speci c clinical features were signi cantly different between high-and low-risk groups. As shown in the alluvial diagram, the majority of patients of GU-like cluster overlapped with patients of high-risk group. The distribution of high and lowrisk groups as well as the difference of risk scores between the two clusters were analyzed (Online Resource 5D, p<0.001, p<2.22 x e -16 , respectively), which indicated that the risk score for one sample can be an indicator of the corresponding gene cluster to some extent. For example, one patient belonging to high-risk group with higher risk score was more likely to be in GU-like cluster, implying that the method of risk model based on 9 GILncRNAs can be substituted for the clustering method. The results of PCA and t-SNE analysis revealed that the method of classi cation based on the risk model was e cient and can well distinguish LGG samples in our research, indicating that the risk model was not only a model to predict prognosis but also a powerful method to distinguish patients with different clinical outcomes and therapy response.

Conclusion
Overall, based on 102 genome-instability related lncRNAs (GILncRNAs), we obtained two clusters with distinct prognosis, immunotherapeutic and PARP inhibitor response, demonstrating that the GILncRNAs were correlated with the above characteristics. Furthermore, the risk score calculated based on 9 genomeinstability related lncRNAs in the risk model was an effective biomarker with powerful predictive function.