1. Identification of lncRNAs related to PD-L1
The genetic information, RNA transcriptome data, and related clinical data of glioma patient samples were acquired from the TCGA and CGGA database(clinical baseline information in Table 1).PD-L1 is encoded by CD274, through the correlation analysis of 13895 lncRNAs with CD274 gene expression, we screened 52 lncRNAs related to CD274 gene expression. Univariate cox analysis and multivariate cox analysis were sequentially used to figure out differential lncRNAs, with p<0.05 as the standard, 23 lncRNAs independently related to CD274 (PD-L1) were screened out(Supplementary Table 1). We performed Lasso regression on these lncRNAs to avoid overfitting and finally filtered out 6 lncRNAs (AL355974.3, LINC01271, AC011899.3, MIR4500HG, LINC02594, AL357055.3) as predictors(Fig. 2 a-c). AL355974.3, LINC01271, and AC011899.3 were positively correlated with the CD274 expression level, while MIR4500HG, LINC02594 and AL357055.3 showed the opposite correlation(Fig. 2 d, e).
Table 1 Clinical characteristics of patients with GBM and LGG
Characteristic
|
N
|
High, N=3361
|
Low, N=3361
|
p-value2
|
Age
|
672
|
54(42,63)
|
38(30,48)
|
<0.001
|
Gender
|
672
|
|
|
>0.9
|
FEMALE
|
|
143(43%)
|
142(42%)
|
|
MALE
|
|
193(57%)
|
194(58%)
|
|
Cancer
|
672
|
|
|
<0.001
|
GBM
|
|
143(43%)
|
7(2.1%)
|
|
LGG
|
|
193(57%)
|
329(98%)
|
|
IDH
|
662
|
|
|
<0.001
|
Mutant
|
|
122(37%)
|
312(94%)
|
|
WT
|
|
207(63%)
|
21(6.3%)
|
|
1p/19q
|
668
|
|
|
<0.001
|
codel
|
|
50(15%)
|
121(36%)
|
|
non-codel
|
|
282(85%)
|
215(64%)
|
|
MGMT
|
635
|
|
|
<0.001
|
Methylated
|
|
184(61%)
|
294(88%)
|
|
Unmethylated
|
|
117(39%)
|
40(12%)
|
|
1Median(IQR);n(%)
2Wilcoxon rank sum test;Pearson’s Chi-squared test
2. Risk score model can predict the prognosis of glioma patients
We observed the effects of the 6 PD-L1-related lncRNAs on the prognosis of glioma. Kaplan-Meier survival curve analysis showed AL355974.3 (HR=4.06, p<0.001), LINC01271 (HR=3.51, p<0.001), and AC011899. 3 (HR=2.66, p<0.001) were significantly related to a poor prognosis; MIR4500HG (HR=0.39, p<0.001), LINC02594 (HR=0.48, p<0.001), and AL357055.3 ( HR=0.28, p<0.001) were relevant to a better prognosis(Supplementary Figure 1). The result indicated that the prediction model based on these 6 lncRNAs can better reflect the prognosis of glioma.
The risk score model based on 6 PD-L1-related lncRNAs’ expressions divided patients with glioma into two groups according to the median risk score. The survival analysis of the TCGA data set showed that in all glioma patients, the prognosis of the low-risk score group was much better than which of the high-risk score group (p<0.0001)(Fig. 3 a-c), but from the perspective of tumour pathological classification, among GBM patients with different risk scores, only PFI is conspicuously different(p=0.018),. In contrast, the three prognostic indicators of LGG patients are all significant differences (OS, DSS: p<0.0001, PFI: p=0.00094)(Fig. 3 d-i). The survival analysis in the CGGA validation data set containing 306 patients has similar results. OS of different risk scorers in either all glioma or separate LGG patients has significant differences (OSALL: P <0.0001; OS LGG: P = 0.00024), but not in GBM patients (p=0.37)(Fig. 3 j-l).
To further explore the prognostic efficiency of PD-L1-related lncRNAs, we use consensus cluster analysis to group samples; it is determined that k=2 is the optimal choice in both the training and verification cohorts. Therefore, in the two data sets, we get two subgroups, Cluster1 and Cluster2(Supplementary Fig. 2 a-c), which can be distinguished by principal component analysis (PCA)(Supplementary Fig. 2 d, e), proving the rationality of clustering. The risk score of the subgroups generated by cluster analysis indicated that the risk score of Cluster1 was conspicuously higher than that of Cluster2 (p<2.2e-16)(Supplementary Fig. 2 f). Survival analysis of cluster subgroups also suggested a poor prognosis for Cluster 1(Supplementary Fig. 2 g-r).
These results indicate that the risk score model based on 6 PD-L1-related lncRNAs has good prediction efficiency in the prognosis of glioma patients, especially in the LGG patient group.
3. Risk score is correlated with the clinical, pathological, and genetic characteristics of glioma
We further explored glioma's clinical and pathological characteristics in different risk score groups. Heat maps and box plots show that in the TCGA data set, MGMT unmethylated, 1p/19q non-codel, IDH wild-type, older (>45 years old), higher WHO grade, classical and mesenchymal subtype patients have increasing CD274 expression level, higher risk score, and be classified in group cluster 1. However, the cancer subtype analysis of GBM patients didn’t find a significant difference in MGMT methylation status, while 1p/19q co-deletion status is undifferentiated in the LGG subtype either (Fig. 4 a, Supplementary Fig. 3).
Genomic heterogeneity leads to prognostic differences in tumours, genetic target-based precision medicine is also one of the goals of PPPM[27]. Higher tumour mutation burden (TMB) in patients with a high-risk score was observed (p<0.001, Fig. 4 b). In the TCGA set, 292 (91.82%) of the 318 patients with high-risk scores and 303 (96.19%) of the 315 patients in another group had somatic mutations(Fig. 4 c, d). Genome variation analysis shows that among the low-risk patients, IDH1, ATRX, TP53, NOTCH1, IDH2, and CIC mutation frequency was higher than that of high-risk patients (p <0.001), while EGFR, PTEN, TTN, RB1, and NF1 gene mutation frequency is higher in a high-risk group (p<0.001)(Fig. 4 e).
4. Risk score correlates with the immune inflammatory microenvironment and immunotherapy responsiveness of glioma
Immune cells and inflammatory cells are an essential part of tumour tissue, the tumour immune microenvironment(TIME)can affect tumour progression and response to treatment. We use GSVA to identify the immune infiltration and inflammation profiles associated with characteristic lncRNAs and draw correlation maps and heat maps.
To clarify whether PD-L1-related lncRNAs affect the prognosis by affecting tumour immunity and inflammation, we compared the degree of immune cell infiltration. In TCGA, the results demonstrated that the macrophages (M0, M1) were conspicuously increased (p<0.001) in the high-risk score group. In contrast, patients in the low-risk score group had significantly increased monocytes (P<0.001) (Fig. 5 a, b, Supplementary Fig. 5 a). Consistently, in the CGGA set, macrophages had a higher abundance and relevance of risk score in the high-risk score group and Cluster 1(Supplementary Fig. 5 b, c).
In terms of inflammation-related molecules, in the TCGA set, the high-risk group and Cluster 1 were associated with higher MHC-I, MHC-II, LCK, HCK, and interferon, especially STAT1, while IgG had increased expression in the low-risk group and Cluster 2(Fig. 5 c, d). The results of the verification show slightly different relevance degrees in 5 molecules(Supplementary Fig. 4 d, e). This suggests that PD-L1-related lncRNAs regulate lymphocyte activation, activation of antigen-presenting cells and interferon signalling in gliomas.
The growth and progression of cancer are in connection with immunosuppression. ICIs have been considered revolutionary immunotherapy for a variety of tumours, including glioma; the tumour microenvironment influences the response of the tumour to immune checkpoint inhibitor therapies. We analyzed the expression of five common immune checkpoints (CTLA-4, HAVCR2, IDO1, LAG3, PDCD1) in patients, and found that people with high-risk scores had higher expression of these five immune checkpoints (P<0.001, Fig. 5 e). We further forecasted the patient's response to anti-PD-1 and anti-CTLA-4 therapies. We figured out that patients in the high-risk group more tend to be sensitive to anti-PD-1 treatment, and the response predicted to anti-CTAL-4 treatments of the two groups is similar(Fig. 5 f). Risk scores based on PD-L1-related lncRNAs can predict response to immunotherapy in glioma patients and thus guide personalized medicine.
For a profound understanding of the possible biological mechanisms of PD-L1-related lncRNAs, the GSEA analysis was applied, and the results showed that PD-L1-related lncRNAs were enriched in biological processes such as neutrophil activation involved in immune response, neutrophil degranulation, and T cell activation; Cell components were mainly concentrated in cell leading edge, the outer side of the plasma membrane and synaptic membrane; Molecular functions were focused on passive transmembrane transporter activity, channel activity and phospholipid binding(Fig. 5 g).
5. Clinical models based on PD-L1-related lncRNAs and clinical characteristics can predict the prognosis of glioma patients
According to the results above, we constructed a clinical model combining risk scores and clinical case characteristics according to the TCGA data set data and made a nomogram. The 3-year and 5-year OS predicted by the model are similar to the actual observations, suggesting that the model has a specific predictive value. Taking account of the nomogram, the patients were divided into two new risk groups. Survival analysis showed a conspicuously different prognosis (P<0.0001) in the two groups. The 3-year AUC under the ROC curve is 0.915, and the 5-year AUC is 0.881. The high specificity and sensitivity prove that compared with the risk score alone, the model that combines clinicopathological characteristics has higher predictive power(Fig. 6).
6. LINC01271 is significantly associated with features of glioma
Previous research on the role of PD-L1-related LncRNA (AL355974.3, AC011899.3, MIR4500HG, AL357055.3, and LINC02594) in tumours are few, and most of them are proteins encoded by their homologous genes. Diermeier's team found that lncRNA Mammary Tumor-associated RNA 25 (MaTAR25) contributed to the malignancy of breast tumour cells in vitro and in vivo. LINC01271 was identified as a homologous human lncRNA of MaTAR25, and they proved that increased LINC01271 expression was relevant to poor prognosis and metastasis[28, 29]. The role of LINC01271 in gliomas has not been reported, so we chose LINC01271 as the target to discover its influence on the development and prognosis of glioma.
The glioma patients in two data sets were separated into two different groups according to LINC01271 expression. In the TCGA set, the overall OS, DSS, and PFI of patients with high expression were worse than the low expression group (p<0.0001), but from the perspective of cancer type, the prognostic difference was only reflected in LGG (LGG: OS p=0.0042, DSS p =0.0083, PFI p<0.0001; GBM: OS p=0.27, DSS p=0.33, PFI p=0.46), which is in line with the results in the CGGA set(Fig. 7 a-c, Supplementary Fig. 5 a-i).
Then we studied the correlation between LINC01271 expression and the characteristics of immunity, inflammation, and clinicopathology. The results showed that the abundance of M1 subtype macrophages was positively correlated with the LINC01271 expression(p<0.01). However, the M0 and M2 subtypes did not show significant differences. (Fig. 7 d, e, Supplementary Fig. 5 j). The same analysis was performed on CGGA database samples, and the heatmap showed similar results(Supplementary Fig. 5 k).
Inflammatory characteristics analysis shows that STAT1, MHC-II, LCK, MHC-I, interferon, and HCK are associated with high LINC01271 expression. In contrast, patients with low LINC01271 expression have higher IgG expression(Fig. 7 f, g). This is consistent with the previous analysis results based on the 6 PD-L1-related lncRNAs.
As for clinicopathological characteristics, we found that CL, ME subtype, MGMT non-methylation, 1p19q, noncodel, IDH wild type, higher age, higher WHO grade, and higher GBM patients proportion are associated with high LINC01271 expression, which illustrates that LINC01271 is a risk factor for prognosis(Fig. 7 d, f).
GSEA analysis revealed that LINC01271 mainly participated in the same biological process as the 6 PD-L1-related lncRNAs model, but concentrated on different cellular components, including the cell-substrate junction, focal adhesion, external side of the plasma membrane and others, and were related to molecular functions like acting binding and phosphoric ester hydrolase activity(Fig. 7 h). The expression of LINC01271 is also in connection with some hallmark gene sets, such as ALLOGRAFT REJECTION, APOPTOSIS, COAGULATION, COMPLEMENT, and E2F TARGETS(Fig. 7 i), suggesting the immune pathways related to characteristic genes, and perhaps related to the mechanism of LINC01271's biological role in glioma.
Finally, we predicted the sensitivity of patients to TMZ treatment, which is currently one of the most common chemotherapy drugs for glioma. With data derived from the GDSC database, we built a predictive model and estimated the half-maximal inhibitory concentration (IC50) value of TMZ in each patient. Unfortunately, there was no statistical difference in IC50 when grouping according to LINC01271 expression (Supplementary Fig. 5 l).
7. Knockout of LNC01271 effectively suppressed the proliferation and migration of glioma cells
We silence LINC01271 gene expression in U251 and U87 to detect proliferation and migration. Rt-qPCR confirmed the silencing efficiency of siRNA, si-LINC01271-334 and si-LINC01271-1564 significantly reduced the expression of LINC01271, si-LINC01271-1196 also had the same trend(Fig. 8 a). CCK-8 assay, plate cloning assay, and Edu assay showed inhibited proliferation of tumour cells and limited colony formation(Fig. 8 b-d). Transwell migration assay suggested inhibition of tumour cell migration after LINC01271 knockout(Fig. 8 e). In conclusion, knockout of LINC01271 visibly suppressed the proliferation and migration of glioma cells, proving that LINC01271 positively regulates the progression of glioma and lncRNA LINC01271 is expected to become a new target for the treatment.