A Prognostic lncRNA Signature Associated with Lipid Metabolism in Glioma


 Background: Long noncoding RNAs (lncRNAs) are promising cancer biomarkers and therapeutic targets . And lipid metabolism reprogramming is a trait of cancer metabolism. However, the relationship between lncRNAs and lipid metabolism and their prognostic value are still unclear in glioma. Hence, we screened for prognostic lncRNAs associated with lipid metabolism and explored their potential biological functions and effects on glioma. Methods: The glioma data were obtained from TCGA and CGGA databases. Correlation analysis was used to identify the lncRNAs associated with lipid metabolism . The Cox and LASSO regression analysis were adopted to find prognostic lncRNAs, and further established ceRNA network, clustering and risk score models. CNV and somatic mutation data were used to explore the correlation between risk score and genomic alterations. GSVA and GSEA was adopted to reveal the potential biological functions of prognostic lncRNAs. The abundance of 22 immune cells was inferred by CIBERSORT algorithm . The TIDE algorithm and GDSC database were used to predict the patient's response to immunotherapy and chemotherapy, respectively. Human glioma cell lines were used for further experimental validation. CCK-8 and colony forming assays were adopted to evaluate the cell viability. The cell proliferation was detected by EdU assay. Transwell assay was used to evaluate the ability of cell invasion and migration. Results: A total of twenty prognostic lncRNAs associated with lipid metabolism were found, and a prognostic lncRNA signature was established. A high risk-score suggested a poor prognosis , more malignant clinicopathological and genomic aberrations features. The risk score was also an independent prognos tic factor. The high-risk patients were more likely to benefit from anti-PD1 treatment. The biological function of signature lncRNAs was mainly to regulate the biosynthesis and transportation of lipid (especially the fatty acid). Among signature lncRNAs, LINC01614 was associated with prognosis, genomic aberrations characteristics, m6A methylation modification, immune infiltration status, and chemotherapy response of patients with glioma. In vitro experiments, silencing the expression of LINC01614 could inhibit the viability, migration , and proliferation of glioma cells. Conclusions: A prognostic lncRNA signature containing twenty lncRNAs associated with lipid metabolism was established in glioma. The signature lncRNAs play an important role in the development of glioma by regulating lipid biosynthesis and transportation. Our study also provides a new perspective for understanding the lipid metabolism and the biological role of lncRNAs in glioma.

detected by EdU assay. Transwell assay was used to evaluate the ability of cell invasion 48 and migration. 49 1 Background 68 As a common primary brain tumor, glioma accounts for 75% of all malignant brain 69 tumors [1,2]. Despite there is a great advance in treatment, the prognosis of patients 70 with glioma remains poor. The 5-year survival rate of glioblastoma (GBM) is only 71 about 5% [3]. Although low-grade glioma (LGG) has a relatively better prognosis, more 72 than half of LGG will progress to high-grade aggressive glioma within 10 years [4]. 73 Therefore, it is important to further explore the mechanism of glioma development and 74 find new therapeutic targets. 75 Metabolic reprogramming is currently considered a marker of cancer. Cancer cells 76 change their metabolism pattern to meet the need for rapid proliferation [5]. This 77 metabolic recombination involves various key pathways of metabolism, such as 78 glucose metabolism (represented by aerobic glycolysis, also known as the Warburg 79 effect), amino acid metabolism, lipid and nucleotide synthesis [6]. In glioma, lipid 80 metabolism abnormality is an important aspect of metabolism reprogramming. The key 81 genes in the regulation of lipid metabolism have been demonstrated to be overexpressed 82 in glioma [7,8], and oncogenic signaling pathways have been revealed to be involved 83 in lipid metabolism reprogramming [7,9]. Besides, glioma cells can de novo synthesize 84 non-essential fatty acids using acetyl-CoA [10]. Compared to the normal brain, the ratio 85 of x-6/x-3 fatty acid is increased in glioma [11]. 86 Long non-coding RNAs (lncRNAs) have no ability to encode proteins, but they are 87 involved in many cellular processes, such as cell development and differentiation, 88 metagenes representing different types of inflammation response [21], and further 145 correlation analysis was conducted. 146

Immunotherapy and chemotherapy response prediction 147
The patients' response to immune checkpoint blocking therapy was predicted by the 148 TIDE [22] algorithm. The patients' sensitivity to temozolomide (TMZ) and gemcitabine 149 therapy was predicted using the data derived from the GDSC database 150 (https://www.cancerrxgene.org). The "pRRophetic" R package was used for the 151 estimation of the half-maximal inhibitory concentration (IC50), which represents the 152 drug response. 153

Transwell migration assay 180
The cells were digested into single cells with pancreatin and suspended in serum-free 181 medium. Then 100ul cell suspension and 500ul of complete medium with 10% FBS 182 were added to the upper and lower chambers of transwell, respectively. After incubation 183 for 48 hours, the upper chamber was washed with PBS and the cells were wiped off. 184 After fixation for 20 minutes with 4% paraformaldehyde, the membrane was stained 185 for 5 minutes with 0.1% crystal violet. After washing with PBS, the membrane was 186 placed on the glass slide and observed under an inverted microscope. The 187 decolorization was made using 10% acetic acid. At 550nm, the OD value was 188 determined by a microplate reader. 189

Colony forming assay 190
The cells were resuspended and seeded in 6-well plates. After 2-3 weeks of culture, 191 the culture solution was discarded and 1 ml of 4% paraformaldehyde was added to each 192 well for fixation for 15 minutes. The fixing solution was removed and 1ml of crystal 193 violet solution was added for staining at room temperature for 30min. After the staining 194 solution was washed off, the colonies were photographed. After decolorization, the OD 195 value at 550nm was measured. 196

Cell proliferation assay 197
The 5-ethynyl-2'-deoxyuridine (EdU) assay kit (RiboBio, China) was used for the 198 evaluation of cell proliferation. Cells were incubated with 100ul of 50μM EdU medium 199 overnight and then fixed with 4% paraformaldehyde. Subsequently, 100ul of 1 × 200 Apollo® staining reaction solution was added for incubation for 30 minutes. After 201 washing with 0.5% TritonX-100, cells were incubated with 100ul of 1 × Hoechst 33342 202 reaction solution for 30 minutes. Finally, cells were observed and photographed with 203 confocal microscopy. And the cell proliferation rate was calculated. 204

Statistical analysis 205
Statistical analysis was conducted by R software (version 3.5.3). Statistical 206 significance between groups was determined using the Student's t-test and two-way 207 analysis of variance (ANOVA). Cox regression analysis, time-related Cox model 208 variables were tested using the proportional hazard hypothesis. Spearman rank was 209 adopted in the correlation analysis. The correlogram and heatmap were drawn using 210 the "corrgram" and "pheatmap" R packages, respectively. The "survival" R package 211 was used for survival analysis. The Kaplan-Meier method was adopted to draw 212 survival curve. The AUC value and ROC curve were obtained by the "survival ROC" 213 package. The R package "rms" was used to establish the nomogram. In this study, p < 214 0.05 was considered statistically significant. 215

A set of lipid-metabolic lncRNAs with prognostic value was identified 217
The workflow designed for the study was shown in Figure 1 shown in Additional file 1: Figure S1. Consequently, a prognostic lncRNA signature 228 associated with lipid metabolism was established. We also build a competing 229 endogenous RNA (ceRNA) network based on these signature lncRNAs ( Figure 2c). 230

Consensus clustering analysis based on signature lncRNAs identified two 231
subgroups with distinct outcomes 232 The samples of glioma patients were grouped using consensus clustering analysis. 233 In the TCGA dataset, we found that k = 2 was an optimal selection, which maintained 234 the maximal consensus within clusters and had a small incremental change in area 235 under the CDF curve (Additional file 1: Figure S2a-S2c). Therefore, two glioma 236 subgroups: cluster 1 and cluster 2, were obtained. Further principal component 237 analysis (PCA) could distinguish the two groups, which proved the rationality of the 238 clustering (Additional file 1: Figure S2d). The prognosis of patients in the two clusters 239 was significantly different. The overall survival (OS), disease specific survival (DSS), 240 and progression free interval (PFI) of patients in cluster1 were all worse than those in 241 cluster 2 (Additional file 1: Figure S2e). The results of consensus clustering also 242 were confirmed in the CGGA dataset (Additional file 1: Figure S3). The above results 243 preliminarily prove the prognostic value of signature lncRNAs. 244

Risk score model had good prognostic predictive ability 245
Based on the signature lncRNAs, we further established a risk score model and 246 calculated the risk score of each glioma patient. According to the median value of the 247 risk score of all patients, we divided the patients into high-risk group and low-risk group. 248 We found that patients in the two groups had distinct patterns of prognosis. Compared 249 with low-risk patients, the high-risk patients had worse DSS in the TCGA pan-glioma, 250 LGG, and GBM cohorts (p < 0.05, Figure 2d). And the poor OS and PFI were also 251 found in the high-risk patients (p < 0.05, Additional file 1: Figure S4a-S4b). The same 252 results were observed in the CCGA dataset (p < 0.05, Additional file 1: Figure S4c

Risk score correlated with clinicopathological features 258
In the TCGA dataset, we conducted a subgroup analysis to explore the correlation 259 between risk score and clinicopathological features. We found that the risk scores of 260 patients within subgroups separated by age, IDH status, LGG and GBM, WHO grade, 261 1p19q codel status, MGMT promoter status, and subtype, were significantly different 262 (p < 0.001, Additional file 1: Figure S5a-S5g). However, no difference was found in 263 gender subgroup (p > 0.05, Additional file 1: Figure S5h). The higher level of risk score 264 preferred to distribute in age ≥ 45, GBM, higher grade, IDH-wt, MGMT promotor 265 unmethylated, 1p/19q noncodel, or classical and mesenchymal subtype patients. 266 Similar results also been found in the GBM and LGG cohort (Additional file 2: Figure  267 S6). These findings suggest that the risk model can help us infer the clinicopathological 268 characteristics of glioma patients. 269

Functional analysis of signature lncRNAs 291
We adopted the GSVA analysis based on GO to further investigate the biological 292 functions of signature lncRNAs. The 20 biological processes of lipid metabolism most 293 significantly correlated with risk score were identified by correlation analysis in the 294 TCGA and CGGA datasets, respectively (Figure 3g, 3i). These lipid metabolism 295 functions were significantly active in high-risk patients as shown in heatmaps ( Figure  296 3f, 3h), and mainly focused on the biosynthesis and transportation of the lipid 297 (especially the fatty acid). In the TCGA dataset, positive regulation of unsaturated 298 fatty acid biosynthetic process, very low-density lipoprotein particle assembly, and 299 negative regulation of lipid biosynthetic process were the biological processes of lipid 300 metabolism most relevant with risk score (Figure 3f-3g). And very low-density 301 lipoprotein particle assembly, long chain fatty acid import across plasma membrane, 302 and positive regulation of intracellular lipid transport were the biological processes of 303 lipid metabolism most correlated with risk score in CGGA dataset (Figure 3h-3i). 304 3.7 The risk score was an independent prognostic factor 305 According to Cox regression analysis, we comprehensively analyzed the effect of 306 risk score and clinicopathological characteristics on clinical outcomes. In the TCGA 307 dataset, we found that age, risk score, and WHO grade were independent predictive 308 factors for OS. Age, WHO grade, risk score, and subtype were independent factors in 309 predicting DSS. And IDH and risk score were independent predictors of PFI. In the 310 CGGA dataset, 1p19q, risk score, IDH, and WHO grade were independent prognostic 311 factors for OS (Additional file 1: Table S1). The above results indicate that risk score 312 is an independent predictor no matter which one is used as the clinical endpoint (OS, 313 DSS, and PFI). This also shows that using the risk score for prognostic judgment is 314 reliable. 315 By integrating risk score and clinicopathological features, we further constructed a 316 nomogram (Figure 4a). We found that the predicted 3-year and 5-year survival rates 317 were close with the actually observed ones in TCGA and CGGA datasets (Additional 318 file 2: Figure S7b). According to the median value of patients' points calculated by 319 nomogram, patients were also divided into two groups, high-and low-risk. In the 320 TCGA dataset, the OS of high-risk patients was worse (p < 0.0001, Figure 4b), and the 321 result was also verified in the CGGA dataset (p < 0.0001, Figure 4d). In the TCGA 322 dataset, the AUC values of the ROC curves were 0.906 at the predicted 3-year OS and 323 0.876 at 5-year OS (Figure 4c). In the CGGA dataset, the AUC values for the 3-and 5-324 year OS were 0.867 and 0.898, respectively (Figure 4e). The AUC values of the 325 nomogram are better than that of the risk model. This demonstrates the nomogram has 326 better predictive power than the risk score model. 327

Risk score was a predictor of immunotherapy response 328
The immune checkpoint pathways are the major mechanism of tumor immune 329 resistance [23], and immune checkpoint blockade has emerged as a promising approach 330 to tumor immunotherapy [23,24]. In this study, we first compared the differences in 331 the expression of immune checkpoint molecules between the high-and low-risk groups, 332 and found that CD724 (PD-L1), CTL4, HAVCR2 (TIM-3), IDO1, LAG3, PDCD1 333 (PD1), PDCD1LG2 (PD-L2), and SIGLEC15 had higher expression levels in high-risk 334 patients (p < 0.001, Figure 4f.). We further predicted the patient's response to anti-335 CTLA4 and anti-PD1 therapy. The results suggested that high-risk patients were more 336 likely to respond to anti-PD1 therapy than low-risk patients (p = 0.02, Figure 4g). But 337 no significant difference was observed in response to anti-CTLA4 therapy between the 338 two groups (p > 0.05, Figure 4g). Moreover, a higher tumor mutation burden (TMB) 339 was observed in high-risk patients (p < 0.001, Figure 4h). 340 3.9 LINC01614 correlated with prognosis, genomic alteration, and m6A 341 methylation of glioma 342 Among signature lncRNAs, LINC01614 has been found to affect the prognosis of 343 cancer [25,26], but its effect on glioma has rarely been reported. Hence, LINC01614 344 was chosen for the next analysis. Based on the median value of LINC01614 expression, 345 we divided the patients into two groups in the TCGA dataset, the high-and low-346 expression groups. Compared with patients in the high-expression group, the outcomes 347 (OS, DSS, and PFI) of patients in the low-expression group were significantly better (p 348 < 0.0001, Figure 5a). This also further confirms the value of signature lncRNAs in 349 judging the prognosis of glioma. Further gene set enrichment analysis (GSEA) [27] 350 revealed the biological functions that LINC01614 may be involved in regulating 351 (Additional file 2: Figure S8a). 352 We also found that the high-and low-expression groups had distinct genomic 353 alterations profiles. The mutations in PTEN (18%), EGFR (17%), SPTA1 (6%), KEL 354 (6%), and PKHD1 (6%) were identified in the high-LINC01614 group, while IDH1 355 (81%), TP53 (49%), and ATRX (34%) mutations were detected in the low-LINC01614 356 group (Additional file 2: Figure S8 b-c). 357 As the most common RNA modifications, N6-methyladenosine (m6A) methylation 358 modification can influence cancer progression by regulating cancer-related biological 359 functions, and noncoding RNAs play an important role in the regulation of these m6A 360 modifications [28]. Therefore, we investigated the relationship between LINC01614 361 and m6A-related genes. The results showed that the expression of ALKBH5, CBLL1, 362 ELAVL1, FTO, HNRNPA2B1, HNRNPC, LRPPRC, RBM15, RBM15B, WTAP, 363 YTHDC1, YTHDF1, YTHDF2, YTHDF3, and ZC3H13 was significantly different 364 between the high-and low-expression groups of LINC01614 (p < 0.05, Figure 5b). This 365 suggests that LINC01614 may participated in regulating m6A methylation in glioma. 366 And mutations of m6A-related genes were shown in Additional file 2: Figure S8d. 367

LINC01614 correlated with immune infiltration and therapy response 368
The tumor immune microenvironment affects tumor development, metastasis, and 369 outcomes [29]. Hence, we also explored the correlation between LINC01614 and 370 immune activity. We first compared the differences in the abundance of 22 immune 371 cells between the LINC01614 high-and low-expression groups. The results of 372 CIBERSORT suggested that the patients with high LINC01614 expression had a higher 373 percentage of resting memory CD4 + T cells, follicular helper T cells, Tregs, gamma 374 delta T cells, resting NK cells, macrophages (M0, M1, and M2), and neutrophils (p < 375 0.05, Figure 5c). While patients in the low-expression group had higher percentages of 376 memory B cells, naive CD4 + T cells, activated NK cells, monocytes, activated dendritic 377 cells, activated mast cells, and eosinophils (p < 0.05, Figure 5c). The correlation 378 analysis indicated that the infiltration of Th2 cells, macrophages, neutrophils, iDC, aDC, 379

T cells, eosinophils, NK CD56dim cells, Th1 cells, T helper cells, cytotoxic cells, NK 380 cells, mast cells, DC, B cells, and Th17 cells was positively correlated with the 381
expression of LINC01614 (p < 0.05, Figure 5d, 5f). In terms of inflammatory response, 382 MHC-II, STAT1, LCK, MHC-I, HCK, and Interferon were positively correlated with 383 the expression of LINC01614. While IgG, a marker for B cells, was negatively 384 correlated (p < 0.05, Figure 5e, 5g). These results suggest that LINC01614 is involved 385 in the regulation of immune microenvironment in glioma. 386 The differences in the expression of immune checkpoint molecules between the high-387 and low-expression groups of LINC01614 were also compared. We found that CD724 388 (PD-L1), CTL4, HAVCR2 (TIM-3), IDO1, LAG3, PDCD1 (PD1), PDCD1LG2 (PD-389 L2), and SIGLEC15 had higher expression levels in high-LINC01614 patients (p < 390 0.001, Additional file 2: Figure S9a). We further predicted the patient's response to 391 anti-CTLA4 and anti-PD1 therapy, but no significant difference was observed between 392 the two groups (p > 0.05, Additional file 2: Figure S9c). Meanwhile, we also predicted 393 the patient's sensitivity to chemotherapy. The results showed that the estimated IC50 394 value of temozolomide (TMZ) and gemcitabine in low-expression patients were higher 395 compared with high-expression patients (p < 0.001, Additional file 2: Figure S9b). This 396 indicates that patients with high LINC01614 may be more sensitive to TMZ and 397 gemcitabine therapy. 398

LINC01614 affected glioma cell viability, migration, and proliferation 399
To further understand the effects of LINC01614 on glioma, we conducted in vitro 400 experiments using glioma cell lines U251 and U87-MG, and LINC01614 was knocked 401 out by using siRNA. As shown in Figure 6a, compared with the control group and 402 siRNA-NC group, the expression of LINC01614 in LINC01614-siRNA group was 403 significantly decreased. This indicates that the expression of LINC01614 in glioma cell 404 lines was successfully interfered with siRNA. We further conducted CCK-8, colony 405 forming, transwell, and EdU assays to explore the effects of LINC01614 on the activity, 406 proliferation, and migration of glioma cells. The CCK-8 assay indicated that cell 407 viability was inhibited by silencing LINC01614 expression (Figure 6b). In colony 408 forming assay, the formation of glioma cell colonies was significantly reduced when 409 LINC01614 was knockdown by siRNA (Figure 6c prognosis, and a high-risk score suggested a poor prognosis. Furthermore, the risk score 437 was identified as a factor with independent prognostic value. These findings show that 438 this lncRNA signature is an ideal prognostic predictor for glioma patients. Besides, we 439 constructed a nomogram, which has better prognosis prediction ability than using the 440 risk model alone, further improving clinical practicality. In addition, the risk score was 441 found to be associated with the malignant clinicopathological features of glioma. In compared with normal brain tissue and low-grade glioma, its expression is elevated in 481 GBM and related to poor prognosis [11]. The high expression of carnitine 482 palmitoyltransferase IA (CPT1A), which is the modulator of long chain fatty acid 483 transport to mitochondria and fatty acid β-oxidation, is associated with GBM 484 progression [53]. In addition, the accumulation of long-chain fatty acid metabolites is 485 related to the poor survival of malignant glioma patients [32]. We also found that very 486 low-density lipoprotein particle assembly, negative regulation of lipid biosynthetic 487 process, positive regulation of unsaturated fatty acid biosynthetic process, positive 488 regulation of intracellular lipid transport, and long chain fatty acid import across 489 plasma membrane may be "onco-biological process of lipid metabolism" in glioma. 490 Immune checkpoint blockade is a promising treatment for glioma patients. However, 491 not all patients can benefit from it [24,54]. Therefore, we wonder whether this lncRNA 492 signature can be used as a predictive marker for the clinical response to immune 493 checkpoint blocking therapy in glioma patients. The results showed that high-risk 494 patients were more likely to benefit from anti-PD1 treatment. This may be related to 495 their high expression of PD1/PD-L1. Also, a higher TMB was found in high-risk 496 patients, and high TMB is currently considered to indicate a relatively better response 497 rate to immunotherapy [55]. 498 As mentioned above, LINC01614 is related to the progression and poor prognosis of 499 cancer, but its role in glioma has been poorly studied. Therefore, we selected 500 LINC01614 among signature lncRNAs for further verification analysis. We found a 501 negative correlation between the clinical prognosis of glioma patients and the 502 expression of LINC01614. This further confirmed the prognostic value of the identified 503 signature lncRNAs to some extent. And LINC01614 may be involved in the m6A 504 methylation modification of glioma. Patients with high-and low-expression of 505 LINC01614 had different genomic aberration, immune infiltration, and inflammation 506 profiles. Macrophage infiltration level, which is positively correlated with the 507 malignancy grade [56], was higher in high-LINC01614 patients. The high-LINC01614 508 patients had more active inflammatory responses related to T cells and macrophages, 509 while the inflammatory response associated with B cells was depressing. Compared 510 with the low-expression group, the expression of immune checkpoint molecules was 511 higher in patients with high LINC01614, but there was no difference in immunotherapy 512 response between the two groups, which may require further expansion of the sample 513 size for research. In addition, patients with high expression of LINC01614 were more 514 likely to benefit from TMZ and gemcitabine treatment. In vitro experiments, the CCK-515 8 and colony forming assays showed that the viability of glioma cells was suppressed 516 by silencing LINC01614 expression. In EdU assay, the proliferation of glioma cells 517 was also reduced with the inhibition of LINC01614 expression. And silencing 518 LINC01614 inhibited the glioma cell invasion in the transwell migration assay. In 519 summary, down-regulation of LINC01614 has a significant inhibitory effect on glioma. 520 The above results indicate that LINC01614 is an important regulatory factor in glioma 521 and may be involved in a variety of biological processes. 522

523
Taken together, we established a lncRNA signature that can be used to predict 524 prognosis and therapy response in patients with glioma. These signature lncRNAs have 525 an important impact on glioma by regulating lipid biosynthesis and transportation. Our 526 findings contribute to further understand the lipid metabolism and the role of lncRNAs 527 in glioma and provide a reference for new therapeutic targets. However, further studies 528 are needed to clarify the specific mechanism of lncRNAs regulating lipid metabolism. 529 Elsherbiny ME, Emara M, Godbout R: Interaction of brain fatty acid-binding 587 protein with the polyunsaturated fatty acid environment as a potential    LINC01614. c Cell cloning was reduced by interfering with LINC01614 expression in 763 colony-forming assay. *p < 0.05, **p < 0.01, ***p < 0.001. 764 assay was used to evaluate the proliferation of U87-MG cells; Cells were stained with 768 EdU (red) and DAPI (blue). 769