3.1 Landscape of genetic variation of PARPs in glioma

After considering a systematic review of published articles about PARPs. A total of 17 PARPs were collected. The PARPs including 13 poly(ADP-ribose) (PARP1, PARP2, PARP3, PARP4, PARP6, PARP8, PARP9, PARP10, PARP11, PARP12, PARP14, PARP15, and PARP16), 2 tankyrase (TNKS, TNKS2) 1 zinc finger (ZC3HAV1) and 1 TCDD (TIPARP). The process of these PARPs involved in many biological processes, such as mRNA modification, transcription of ribosome protein genes, and protein homeostasis (Fig. 1A). In the glioma genome, the mutation rate for all PARPs is relatively low. Of the 885 (1.92%) samples, a total of 17 samples had genetic changes in PARPs, including missense mutations and splicing site mutations. Among them, the highest mutation frequency is PARP4, the followed are PARP12 and TNKS2 (Fig. 1B). The copy number variation (CNV) of PARP11, ZC3HAV1, and PARP12 is mainly manifested as copy amplification; and TIPARP, PARP8, and TNKS monetize missing copy numbers (Fig. 1C). The location of CNV alteration of PARPs on chromosomes was shown in Fig. 1D. The protein–protein interaction (PPI) network performed by string demonstrates a wide range of interactions between these PARPs (Fig. 1E, S1B). Additionally, the analysis of PARPs expression differences between LGG and GBM revealed significant variations in most PARPs across different grades of gliomas, underscoring their potentially significant role in the malignant progression of gliomas (Fig. 1F). Based on the findings in Fig. 1F, we identified genes displaying significant expression differences between LGG and GBM. In addition, we tested the RNA expression levels of these genes in gliomas of different grades by qRT-PCR assays. Fig. S2 illustrates that 14 PARPs exhibited differential expression patterns across normal brain tissues, LGG, and GBM tissues (Fig. S2A-N). These findings revealed significant differences in epigenetic and gene expression of PARPs between LGG and GBM.

Fig. 1

Genetic Variation Landscape of PARPs in Glioma. A Functional roles of 17 PARPs in cellular biology. B Mutation frequency of 17 PARPs in the TCGA cohort. C Copy number variation (CNV) frequency of PARPs in the TCGA cohort. D Chromosomal location of CNV alterations of PARPs across 23 chromosomes using the TCGA cohort. E Protein–Protein Interaction (PPI) network among PARPs. F Box plot illustrating the expression levels of the 17 PARPs between Lower Grade Glioma (LGG) and Glioblastoma (GBM) in the TCGA cohort. *p < 0.05, ***p < 0.001

3.2 The modification patterns of 17 PARPs in glioma

The two CGGA datasets, CGGA1 and CGGA2 (supplementary table S1), were integrated into a single cohort. Relevant survival data and clinical information were compiled. Prognostic assessment of 17 PARPs in glioma patients was conducted using the univariate Cox regression model (Fig. S1C). Utilizing the expression of the 17 PARPs, PCA analysis distinguished three PARP modification patterns. The results demonstrated the effectiveness of these patterns in distinguishing the samples (Fig. 2A and B). Survival analysis indicated a notable survival advantage associated with the PARP cluster B modification pattern, whereas cluster A exhibited the worst prognosis (Fig. 2C). Using the “ConsensusClusterPlus” R package, patients were categorized into distinct PARP modification patterns based on the expression profiles of 17 PARPs. Unsupervised clustering identified three unique modification patterns: 319 cases in pattern A, 258 in pattern B, and 388 in pattern C (Fig. 2A and Fig. S3A–H). The “ConsensusClusterPlus” R package was utilized to classify patients based on qualitatively different PARP modification patterns derived from the expression of these 17 PARPs. Unsupervised clustering revealed three distinct modification patterns: 319 cases in pattern A, 258 cases in pattern B, and 388 cases in pattern C (Fig. 2A and Fig. S3A–H). These patterns were termed as PARP clusters A–C, respectively (Fig. 2D and supplementary table S3).

Fig. 2

PARPs Modification Patterns and Relevant Biological Pathways. A By consensus clustering, all glioma patients were divided into 3 clusters. B Principal Component Analysis (PCA) was employed to validate the presence of three discernible patterns, determined by the expression profiles of the 17 PARPs across a total of 2228 glioma samples. C Survival curves of glioma samples categorized into the three PARPs modification patterns within the combined glioma cohorts. D Heatmap illustrating consensus clustering based on 17 PARPs. E, F GSVA analysis displaying the enrichment of gene sets in the three patterns

3.3 Characteristics of tumor microenvironment cell infiltration in different modification patterns

We employed GSVA enrichment analysis to discern biological disparities among various PARPs modification patterns. Our findings demonstrated that PARP cluster-A exhibited enrichments in stromal and carcinogenic activation pathways. Conversely, PARP cluster-C displayed a notable association with immunosuppressive biological processes (Fig. 2E, F, and supplementary table S4). Further analysis using ssGSEA and TME revealed that cluster A comprised a higher abundance of natural killer cells, macrophages, eosinophils, and mast cells, indicating a strong correlation with matrix activation. Cluster B, in contrast, exhibited an immune-inflamed phenotype, distinguished by adaptive immune cell infiltration and pronounced immune activation. Meanwhile, PARP cluster C was identified as an immune-desert phenotype, characterized by significant immune suppression (Fig. 3A, B, and S3I). Previous research has elucidated the role of ADP‐ribosylation in the differentiation and functions of both innate and adaptive immune cells [33]. To investigate the role of PARPs within the TME, we performed Spearman's correlation analysis. The results revealed a positive correlation between the expression of genes such as ZC3HAV1, TIPARP, PARP4, PARP14, and PARP15 and the abundance of various immune cell types in gliomas (Fig. 3C). These results emphasize the connection between PARP expression and immune infiltrates.

Fig. 3

Association of PARPs Modification with TME Cell Infiltration. A Abundance of TME infiltrating cells in three PARPs modification patterns. B Differences in stroma-activated pathways, including EMT and TGF-beta pathways, among three distinct PARPs modification patterns. Statistical differences among the three modification patterns were tested using the one-way ANOVA test. C Correlation between each TME infiltration cell type and each PARP using Spearman analyses. Negative correlation marked with blue and positive correlation with red. D Venn diagram showing 5737 PARPs-related genes. E Functional annotation of the PARPs-related genes between three patterns in the combined glioma cohort. *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001

3.4 Generation of PARPs signatures and functional annotation

We identified 5737 DEGs between different PARPs modification patterns using the “limma” package. The Venn diagram shows these results (Fig. 3D, supplementary table S5). GO enrichment analysis showed a significant correlation between biological processes such as PARPs modification and immune system regulation (Fig. 3E, F). To further screen for genes that have an impact on the overall survival of gliomas, we screened 4377 genes of the PARPs phenotype from 5737 DEGs by univariate cox regression and performed consensus clustering analysis based on these genes. Three subgroups were identified and named PARP gene cluster A-C (Fig. S4A-H, supplementary table S6). Unsurprisingly, the expression of PARPs in the 3 PARP gene clusters has obvious differences (Fig. 4A). We also repeated the above studies in the TCGA dataset and identified 3 gene clusters. And the survival analysis results between these three groups are also similar to those in the CGGA dataset (Fig. S5A).

Fig. 4

Construction of PARP Signatures. A Expression levels of 17 PARPs in three gene clusters. B Abundance of each TME infiltrating cell in three gene clusters. C Differences in stroma-activated pathways. D Survival analysis of glioma patients belonging to the three PARP-related gene clusters. E Heatmap illustrating consensus clustering based on PARPS-related differential genes. (F, G) Boxplots demonstrating the PARPscore for PARP patterns (F) and gene clusters (G). H Sankey diagram depicting the association between PARP patterns, gene clusters, PARPscore groups, and fustat. *p < 0.05; **p < 0.01; ***p < 0.001; ***p < 0.0001

3.5 Clinical and relevant pathological features in PARP-associated phenotypes

We conducted ssGSEA enrichment analysis and observed distinct enrichment patterns associated with immune activation, stromal activation, and carcinogenic pathways in PARP gene clusters C, B, and A, respectively. We validated these findings using consistent TME analysis methods (Fig. 4B and C). We also analyzed the differences in survival of these three gene clusters, and the results showed that the cluster C had the best prognosis (Fig. 4D). In addition, we used heatmap to display gene expression patterns between different clinical, and molecular phenotypes groups, and found clear differences in gene expression between each subgroup (Fig. 4E). Given the intricate nature of PARPs modification in gliomas, we introduced a scoring system termed PARPscore to evaluate PARPs modification patterns in glioma patients. PARPscore exhibited significant variations among different PARP gene clusters. Notably, PARP cluster A displayed the highest median score, while PARP clusters B and C showed no noticeable difference, potentially attributed to our dataset and PARP patterns. However, PARP gene cluster A exhibited the highest PARPscore, while cluster C had the lowest (Fig. 4F and G). The alluvial diagram effectively portrayed the attribute changes of individual patients (Fig. 4H). We stratified samples into high and low groups based on the PARPscore, revealing a significant survival advantage for patients with low PARPscore (p < 0.001, Fig. S4I). Univariate Cox regression (Fig. S5C and E) and multivariate Cox regression analyses (Fig. S5D and F), considering factors such as age, grade, IDH status, 1p19q status, etc., indicated that PARPscore is a valuable factor for assessing patient prognosis in both CGGA (Fig. S5C and D) and TCGA datasets (Fig. S5E and F). We further examined the differences in the overall PARPscore proportion among different survival status groups of patients. The analysis revealed higher overall PARPscore in the deceased group, whereas the PARPscore was lower in the surviving group (Fig. 5A, B, C, and D; A and C (CGGA), B and D (TCGA)). Additionally, we investigated the survival rates in patients with different 1p19q status (Codel and Non-codel; Fig. 5E and F), different IDH status (Mutant and Wildtype; Fig. 5G and H), with or without radiotherapy (YES and NO; Fig. 5I and J), and with or without chemotherapy (YES and NO; Fig. 5K and L). Overall, a high PARPscore was associated with the poorest survival rates. In the TCGA dataset, patients with low PARPscore also exhibited a significant survival benefit (p < 0.001, Fig. S5B).

Fig. 5

Relationship of PARP-modified Features to Tumor Somatic Mutations. A–D Differences in the proportion of overall PARPscore in different survival status groups of patients in CGGA (A, C) and TCGA (B, D). E–L Survival rates of patients with different 1p19q status (Codel and Non-codel; E and F), different IDH status (Mutant and Wildtype; G and H), with Radiotherapy or not (YES and NO; I and J), and with chemotherapy or not (YES and NO; K and L). M Differences in PARPscore between patients with high and low TMB. N Scatter plot depicting the relationship between PARPscore and TMB in glioma samples (R = 0.54; p < 2.2e−16). O The survival analyses of glioma samples in high and low TMB groups (p < 0.001). P–Q Waterfall plot presenting tumor somatic mutations in low (P) and high (Q) PARPscore subgroups

3.6 Relationship of PARPs-modified features to tumor somatic mutations

We utilized the "maftools" R package to analyze the distribution differences in somatic mutations between low and high PARPscore groups. Remarkably, we found that tumors with a low PARPscore were linked to reduced TMB (Fig. 5M and N). Notably, the low PARPscore group displayed a broader range of mutations compared to the high PARPscore group (Fig. 5P and Q). Additionally, our analysis indicated that patients with a combination of high PARPscore and high mutation burden experienced the poorest survival rates (Fig. 5O, Fig. S4J). These findings offer a novel perspective for investigating the relationship between PARP modification patterns and somatic mutations in tumors.

3.7 The role of PARPs modification patterns in immunotherapy

CTLA-4 and PD-1 blockade therapies have demonstrated significant improvements in survival rates across various cancer types [34, 35]. Building on the strong association observed between PARPscore and tumor immune response, we analyzed its correlation with TIDE scores. Our results revealed significantly lower TIDE scores in the high-PARPscore subgroup across the CGGA-merge, CGGA1, CGGA2, and TCGA cohorts (Fig. 6A–D). These findings highlight the substantial impact of PARPs on the tumor immune response in gliomas. To further examine the clinical relevance, we assessed the relationship between PARPscore and responses to anti-PD-L1/PD-1 immunotherapy in the IMvigor210 and GSE79671 cohorts. Patients with high PARPscore demonstrated poorer outcomes, with notably shorter survival times (GSE79671, p = 0.001; IMvigor210, p < 0.001, Fig. 6E and I). Patients in the high-PARPscore subgroup showed reduced responsiveness to PD-1/PD-L1 immunotherapy (Fig. 6F, G, J, and K). Furthermore, this subgroup exhibited significantly lower PD-L1 expression in both immunotherapy cohorts, providing insights into the potential molecular mechanisms influencing clinical immunotherapy responses (P < 0.05; Fig. 6H and L). The results demonstrate a significant association between PARP modification patterns and the immune response in gliomas. Furthermore, the PARPs signature accurately predicts clinical responses to anti-PD-1/PD-L1 immunotherapy.

Fig. 6

The Impact of PARP Modification Patterns on Immunotherapy. A–D Distribution of TIDE scores between high and low PARPscore subgroups in the gathered glioma cohorts (p < 2.22e−16) (A), as well as CGGA1 (p < 2.22e−16) (B), CGGA2 (p < 2.22e−16) (C), and TCGA (p = 0.00032) (D) datasets. E, I Survival analyses of OS for the low and high PARPscore subgroups in the two anti-PD-L1 immunotherapy cohorts (GSE79671, p = 0.001 (E) and IMvigor210, p < 0.001 (I)). F, J Proportion of patients responding to PD-1 blockade immunotherapy in the low and high PARPscore subgroups of the GSE79671 (F) and IMvigor210 (J) cohorts. G, K Clinical response to anti-PD-1 immunotherapy based on high or low PARPscore of patients in the GSE79671 (G, p = 0.0086) and IMvigor210 (K, p < 1.5e−07) cohorts. H, L Differences in PD-L1 expression between the PARPscore subgroups in the GSE79671 (H, p = 0.02) and IMvigor210 (L, p < 6.2e−05) cohorts