Total RNA was extracted from TMZ-sensitive/resistant U87MG cell lines by TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. RNA quality was assessed by Nanodrop2000 and Qubit 3.0. RNA integrity was determined by Agilent 2100 Bioanalyzer. mRNA Capture Beads (Vazyme Biotech, Nanjing, China) were used to eliminate rRNAs and a VAHTS Total RNA-Seq Library Preparation Kit (Vazyme Biotech) was used to prepare libraries. Sequencing was performed on the Illumina Hiseq 2500 platform (pair-end 150 bp). The raw sequencing data have been deposited in the NCBI (BioProject accession: PRJNA768121). Data were further processed by R (version 4.1.0). A heatmap of gene expression profiles was generated using the pheatmap package. Principal component analysis (PCA) of each sample was performed and the top two principal components were shown. Differentially expressed genes (DEGs) were identified using the limma package (| log2(fold change) | > 1 and adjusted p-value < 0.05 as threshold). A volcano plot was illustrated to visualize the distribution of DDR genes using the EnhancedVolcano package.

Two independent datasets of glioma patients were collected from the publicly available The Cancer Genome Atlas (TCGA) (https://xenabrowser.net/datapages) and Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn) (13). TCGA and CGGA cohorts enrolled 557 and 656 glioma specimens, respectively. Patients lacking complete clinical annotations were excluded from subsequent analyses. Gene expression data from both cohorts were log2(FPKM+1) transformed.

Univariate cox regression analysis was performed to screen the DDR-related DEGs between TMZ sensitive and resistant glioma cells. DEGs with significance in univariate cox regression were subsequently analyzed by LASSO regression to further select variables. To improve the accuracy of the risk model, LASSO regression was implemented with 10-fold cross validation and run for 1000 rounds to alleviate overfitting effects. Finally, nine DDR genes were picked out to construct the DDR score following the formula: DDR score = ∑ 1 n Coef  ( i ) ∗ exp ( i ) , with Coef (i) meaning the coefficients of each variable and exp (i) representing the expression level of genes.

Glioma patients were divided into high and low DDR score groups according to the cutoff value. Scatter diagrams and Kaplan-Meier curves were generated by the survival package to evaluate the difference in clinical outcome between two DDR score phenotypes. ROC curves were depicted to determine the predictive capability of DDR score using the timeROC package. DDR scores of different glioma grades and pathological statuses were displayed using the ggplot2 package. Univariate and multivariate cox analyses were performed to investigate the prognostic value of DDR score and other clinicopathological indicators. The forest plot delineated the corresponding hazard ratio (HR) and p-value. A nomogram was constructed based on DDR score and other independent prognostic indicators to predict the 1, 2, and 3 years survival probability of glioma patients. Calibration curves were produced to evaluate the utility of the nomogram.

DEGs between high and low DDR score groups were identified using the limma package (| log2(fold change) | > 1 and adjusted p-value < 0.05 as threshold). Symbols of DEGs were extracted using the clusterProfiler package to explore Gene Ontology (GO) and KEGG terms enrichment (14). Gene set enrichment analysis (GSEA) was used to investigate the enriched gene sets based on the fold changes of all genes. The significances were ranked by normalized enrichment score (NES) and adjusted p-value. Gene set variation analysis (GSVA) is a non-parametric, unsupervised approach to estimating the variation of gene set enrichment based on expression profiling. KEGG and HALLMARK gene sets were downloaded from the Molecular Signatures Database (MSigDB) and signatures of tumor metabolism were collected from a published study (15). The GSVA package was exploited to quantify the pathway activity by calculating the GSVA score of each sample (16).

Tumor purity is diluted by the non-tumor components in the TME including immune infiltrates, stromal cells, blood vessels, and extracellular matrix. The ESTIMATE package was applied to calculate the immune score, stromal score, and tumor purity of glioma specimens (17). CIBERSORT was used to evaluate the proportions of 22 immune cell types based on deconvolution methods (18). ssGSEA was employed to investigate the levels of 28 immune cell types based on the marker gene signature score. TIMER was utilized to estimate the percentages of six immune cell types by linear least square regression (19). MCP-counter was implemented to quantify the absolute abundance of eight immune and two stromal cell populations based on the mean level of marker gene expression (20). Gene expression data with standard symbol annotation were imputed to the algorithms above for further analyses.

Glioma specimens were obtained from patients who underwent surgical resection in Shanghai East Hospital (from 2019 to 2021). All participants signed written informed consent for molecular studies before sample collection. The clinical data of patients were recorded with the approval of the Ethical Committee and Institutional Review Board of Shanghai East Hospital.

Eight glioma samples were fixed by immersion in 10% formalin solution and then embedded in paraffin. 10-µm thick tissue sections were deparaffinized and rehydrated (xylene × 2 for 10 minutes each, 100%, 95%, and 75% ethanol for 5 minutes each and deionized water for 5 minutes). The sections were incubated with 3% hydrogen peroxide in methanol for 10 minutes to quench peroxidase activity. Antigen retrieval was performed by boiling sections in 10 mM sodium citrate buffer (pH 6.0) for 10 minutes. After being rinsed with PBS, sections were blocked with normal goat serum for 20 minutes. The samples were incubated with primary anti-PD-L1 (1:200, ab237726, Abcam, Cambridge, UK) or anti-TGFβ antibody (1:200, BA0290, BOSTER, Wuhan, China) overnight at 4°C. The sections were then incubated with secondary antibody () for 30 minutes at room temperature. The staining was developed using 3,3’-diaminobenzidine (DAB) as substrate and counterstained with hematoxylin. The sections were developed using 3,3’-diaminobenzidine (DAB) as substrate and counter-stained with Mayer’s Hematoxylin. All sections were independently reviewed by two pathologists according to the WHO criteria.

Somatic mutation files (SNPs and small INDELs) of TCGA glioma were downloaded from the UCSC Xena browser. The maftools package was used to present the mutational patterns of glioma specimens by the oncoplot function. The mutual exclusivity and co-occurrence of top frequent mutations were delineated by the somaticInteractions function. The mutation load of each specimen was calculated by the tmb function.

Genomics of Drug Sensitivity in Cancer (GDSC) is a public dataset containing information on drug sensitivity in cancer cells and molecular markers of drug response (21). Using the oncoPredict package, GDSC2 gene expression profile and corresponding drug response information were downloaded to generate a ridge regression model that can be applied to glioma transcriptomic data (22). Then the sensitivity scores were yielded to predict the half-maximal inhibitory concentration (IC50) of all drugs in glioma patients.

IMvigor210 is a cohort of 348 urothelial cancer patients treated with PD-L1 blockade therapy (23). The gene expression profiles, tumor mutation burden, neoantigen information, therapeutic responses, and survival data were downloaded using the IOBR package (24). Tumor-intrinsic signatures were derived from the IOBR package and enrichment scores were calculated by ssGSEA algorithm.

Statistical significance of normally distributed variables between two groups was analyzed by unpaired Student’s t-test and nonnormally distributed variables were examined by Wilcoxon test. The Kaplan–Meier survival curve was graphed to explore survival distributions. Log-rank test was used to determine statistical significance between groups. The Pearson correlation analysis was used to test the association between continuous variables. χ²-test was employed to analyze contingency tables. The multi-omics data were standardized by z-score scaling. All statistical analyses were conducted using R (version 4.1.0) and two-sided p-values less than 0.05 were considered statistically significant and labeled as *, p < 0.05; **, p < 0.01; ***, p < 0.001; and ****, p < 0.0001.

To investigate the underlying mechanisms of TMZ resistance in glioma cells, we used RNA sequencing to analyze the transcriptomic alternations in TMZ-resistant U87-MG (U87-MGR) cells. U87-MG and U87-MGR cells exhibited distinct transcriptomic traits and principal components (Figures 1A, B). KEGG analysis revealed that DDR pathways were differentially enriched between TMZ sensitive and resistant cells including DNA replication, base excision repair (BER), MMR, Fanconi anemia, and cell cycle (Figure 1C). We gathered a list of 608 genes regulating DNA replication and/or DNA repair from Molecular Signatures Database (MSigDB) to screen out the DDR-related DEGs. The overall distribution was shown in a volcano plot and the red dots represented the DEGs (Figure 1D).

To establish a DDR risk model, we selected nine genes from the DDR-related DEGs using univariate cox regression and LASSO regression analyses (Figure 2A and Supplementary Figure 1A). We constructed the DDR score based on the expression levels of nine genes and corresponding coefficients (Figure 2B): DDR score = exp (UNG) × 0.007710872 + exp (GINS4) × 0.123645507 + exp (CHAF1B) × 0.164960931 + exp (S100A11) × 0.199658898 + exp (FANCA) × 0.003032437 + exp (USP43) × -0.058596477 + exp (GADD45G) × -0.011548864 + exp (POLR2F) × -0.157450788 + exp (ERCC5) × -0.189557601. A network plot showed the mutual relationship of these nine genes (Figure 2C). An optimized cut-off DDR score of -0.37 was used for classification into high and low DDR score glioma patients. The scatter plot and Kaplan-Meier curves revealed that a high DDR score was associated with a poorer outcome in the TCGA training cohort (Figures 2D, E). ROC analysis indicated that the area under the curve (AUC) for 1, 2, and 3 years survival were 0.875, 0.907, and 0.910, respectively (Figure 2F). Similar survival analysis results were observed in the TCGA validation cohort (Figures 2G, H) and the corresponding AUC values were 0.873, 0.912, and 0.925 (Figure 2I). Furthermore, the AUC for 1, 2, and 3 years survival were 0.872, 0.906, and 0.912 in the entire TCGA cohort (Supplementary Figure 1B). Then we applied CGGA data for external validation, patients in the high DDR score group also had disadvantageous survival compared with counterparts (Figures 2J, K). The AUC values for 1, 2, and 3 years survival were 0.776, 0.813, and 0.807, respectively (Figure 2L). Taken together, survival analysis underlined the robust value of DDR score for predicting glioma prognosis.

Then we assessed the DDR score of patients with different malignancy grades and pathological statuses in TCGA and CGGA. High grade (grade IV) tumors exhibited higher DDR scores than low grade (grade II and III) ones, and gliomas with IDH wildtype and 1p19q non-deletion showed the highest DDR score compared with other clinicopathological subtypes (Figures 3A–D). Moreover, the Kaplan-Meier curves revealed that high DDR score group showed significantly shorter lifespan and worse outcome in IDH mutant lower-grade gliomas (Supplementary Figures 2A, B). We evaluated the prognostic value of DDR score and clinicopathological parameters in TCGA using univariate cox regression analysis. The HR and p-value of each factor were shown in a forest plot (Figure 3E). Variables with a significant p-value were taken into a multivariate model and the results indicated that DDR score and age were independent prognostic factors in the TCGA cohort (Figure 3F). Univariate and multivariate cox analyses also revealed that DDR score, tumor grade, and 1p19q codeletion status independently predicted the prognosis of patients in the CGGA cohort (Figures 3G, H). In addition, a nomogram was constructed by integrating the DDR score and other parameters to predict the 1, 2, and 3 years survival probability of glioma patients in TCGA (Figure 3I). The corresponding calibration curves were close to the ideal model (Figures 3J–L). The nomogram and calibration curves of the CGGA cohort were present in Supplementary Figures 2C–F.

To understand the underlying biology contributing to the extraordinary predictive ability of the DDR score, we explored the GO enrichment in TCGA using the clusterProfiler package. T cell activation was remarkably connected to a high DDR score while signal transduction-related genes were overexpressed in the low DDR score subset (Figures 4A, B). KEGG and GSEA analyses revealed that pathways involved in immune system disorders were significantly associated with high DDR score and synaptic transmission pathways were enriched in the low DDR score group (Figures 4C–F). Then we comprehensively investigated the HALLMARK and KEGG signature scores of each sample using the ssGSEA algorithm. The significances were ranked by adjusted p-values using the limma algorithm and the top 20 enrichment were visualized in heatmaps (Figures 4G, H). High DDR score gliomas exhibited the enrichment of multiple immune activation pathways including IL-6/JAK/STAT3 pathway, interferon-gamma response, inflammatory response, and antigen processing and presentation. Low DDR score correlated with the enrichment of the Wnt signaling pathway and Hedgehog signaling pathway. In addition, we observed the different metabolic regulations between two phenotypes so we further explored the potential metabolic mechanisms exploiting relevant signatures. The results underlined that cyclooxygenase arachidonic acid metabolism, glutathione metabolism, and pentose phosphate were upregulated in the high DDR score group, while alanine aspartate and glutamate metabolism, fatty acid biosynthesis, and sirtuin nicotinamide metabolism were stimulated in the low DDR score phenotype (Figure 4I).

The tumor is surrounded by a dynamic microenvironment that consists of various types of infiltrating immune cells. The interplay between these immune components and tumor cells can shape tumor immunogenicity and affect the tumor response to checkpoint inhibitors (6). Considering the distinct immune and metabolic characteristics of high and low DDR score subsets, we then investigated the association between DDR score and immune infiltration. The stromal score, immune score, and tumor purity of glioma specimens were inferred using the ESTIMATE package. High DDR score correlated with elevated stromal and immune scores but decreased glioma purity (Figure 5A), indicating the raised levels of stromal and immune cells in the TME of high DDR score glioma. We then further explored the TME landscape of two DDR score subtypes using CIBERSORT, ssGSEA, TIMER, and MCP-counter algorithms (Figures 5B–E). Tumors with high DDR scores correlated with elevated levels of multiple immune infiltration. Interestingly, we discovered that both pro-tumor (myeloid-derived suppressor cells (MDSCs), regulatory T cells (Tregs), M2 macrophages, and immature dendritic cells) and anti-tumor (CD8+ T cells, CD4+ T cells, activated dendritic cells, and M1 macrophages) immune cells were up-regulated. Notably, tumors with the presence of CD4+ and CD8+ T cells were defined as immune-inflamed and usually associated with a positive response to ICB. Activated dendritic cells and M1 macrophages are effective in T cell activation through antigen presentation. On the other hand, MDSCs can weaken the activity of effector T cells, mediate the differentiation of Tregs, and promote an immunosuppressive phenotype in macrophages. Immature dendritic cells induce an immunosuppressive TME through expanding Tregs. Further, M2 macrophages contribute to tumor immune evasion by expressing anti-inflammatory cytokines and attenuating the activity of CD8+ T cell (25). In a nutshell, the complicated TME of glioma was characterized by the mixture of pro- and anti-tumor cells, as well as the coexistence of immune activation and suppression.

Given the differential infiltrating immune cells in two DDR groups, we sought to explore the signatures of TME associated with DDR score. All the defined gene signatures were obtained from previous publications (23, 26–30). DDR-related processes were overexpressed in high DDR score specimens, indicating an overexpressed DDR phenotype (Figure 6A). Antigen processing, chemokines, and interferon responses were significantly up-regulated in the high DDR score group, suggesting the enhanced efficiency for T cells to recognize antigens and the triggered inflammation and antitumor immunity. TMEscore was a novel biomarker with high sensitivity in predicting immunotherapy efficacy, and we identified its positive correlation with DDR score. Intriguingly, high DDR score gliomas also showed promoted TGFβ pathway activity which was associated with immunosuppression (Figures 6A, B). These observations revealed the coexistence of anticancer immunity activation and immune suppression in the glioma microenvironment.

Immune checkpoints are regulators of immunological tolerance that function to protect the cells from indiscriminate attack. The activation of inhibitory checkpoint molecules prevents tumors from damage and attack so they can serve as promising targets for cancer immunotherapy (31). We investigated the immune checkpoint expression in glioma specimens and uncovered that CD274 (PD-L1), PDCD1, PDCD1LG2, CTLA4, CD276, HAVCR2, and LAG3 were significantly overexpressed in the high DDR score subtype (Figure 6B). The cGAS-STING pathway is an important component of the cellular innate immune system that functions by detecting cytosolic DNA fragments and consequently triggers cytosolic immunity. Stimulating cytosolic immunity is a state-of-the-art strategy to optimize ICB therapy efficacy by promoting infiltrating T cells to turn immunologically cold tumors into hot tumors (32). Our analysis suggested that the high DDR score subset presented the elevated levels of cGAS-STING pathway members (Figure 6C). Using immunohistochemistry (IHC) analysis, we evaluated the protein expression of PD-L1 and TGFβ in primary and recurrent gliomas. Recurrent gliomas, usually concomitant with chemotherapy and/or radiotherapy resistance, were considered to be in a DDR proficiency status. IHC staining showed that PD-L1 and TGFβ were up-regulated in recurrent tumors compared to counterparts (Figure 6D). Collectively, we identified a collection of differentially expressed TME gene signatures, which potentially indicated an encouraging sensitivity of ICB in high DDR score gliomas.

Tumor mutation burden (TMB) is defined as the number of genetic mutations in a tumor. High tumor mutation burden usually correlates with a positive response to ICB therapy substantially because increasing mutation load potentially generates neoantigens to enhance tumor antigenicity and immunogenicity (33). To outline the somatic mutation landscape of two DDR subtypes, we displayed the top 20 frequent mutations in both groups, (Figures 7A, B). 18% of high DDR score patients exhibited IDH1 mutation (Figure 7A) but conspicuously the alternation rate of IDH1 was up to 92% in low DDR score phenotype (Figure 7B). Therefore, DDR score seemed to be a robust indicator of IDH1 mutation status in glioma. Further, DDR scores of IDH1 mutant patients were notably lower than that of IDH1 wildtype patients (Figure 7C), which was consistent with our results in Figures 3B, D. Meanwhile, ATRX mutant gliomas also showed significantly decreased DDR score (Figure 7D). Conversely, EGFR and PTEN mutations were crucially correlated with the increase of DDR score (Figures 7E, F). These results suggested that EGFR and PTEN mutations may potentiate DDR activity in glioma.

Additionally, higher frequencies of mutation co-occurrences were observed in high DDR score gliomas (Figures 7G, H). Subsequently, we calculated the TMB of each glioma specimen and identified that the DDR score was positively associated with the mutation burden (Figures 7I, J). Nevertheless, no such correlation was observed in IDH wildtype glioblastoma (GBM) (Supplementary Figure 3). Based on the median TMB and the DDR score cut-off value, we divided patients into four groups and uncovered that patients with high DDR scores and high TMB associated with the worst prognosis and those with low DDR scores and low TMB showed the longest survival (Figure 7K).

Immunogenic cell death prompted by certain chemotherapy agents can be exploited to sensitize tumors to checkpoint blockade, so the optimal combination of chemotherapy and immunotherapy warrants further exploration (34, 35). Since the DDR score was generated based on the DEGs of TMZ-resistant glioma cells, we supposed that chemotherapy status potentially correlated with DDR score levels. Our analysis revealed that DDR scores were higher in patients who underwent chemotherapy (Figures 8A, B). OncoPredict package was utilized to predict the sensitivity scores of drugs in high and low DDR score groups and the sensitivity score was positively correlated with the IC50 value of chemotherapy agents. We compared the estimated TMZ sensitivity between two subtypes and found no significance (Supplementary Figure 4A). In addition, MGMT and BCL3 were up-regulated in the high DDR score group but no statistical significance of ALKBH2 expression was observed (Supplementary Figure 4B). Therefore, more studies were needed to investigate the correlation between DDR score and the susceptibility of alkylating agents. Further analysis suggested that targeting the PI3K pathway (CZC24832 and VSP34_8731) and inducing apoptosis (Entospletinib) may be efficient strategies for high DDR score patients (Figures 8C–E). These predictions were hardly surprising because apoptosis and mTORC1 signaling were enriched in the high DDR score phenotype according to the aforementioned functional analysis. Meanwhile, I-BRD9 (BRD9 inhibitor), BIBR-1532 (telomerase inhibitor), and Linsitinib (IGF-1R inhibitor) were candidate drugs for the treatment of low DDR score tumors (Figures 8F–H).

A growing body of early-phase clinical trials have been developed to evaluate the potential of combining DDR-targeted therapy with ICB (9). Currently, most reports regarding DDR-immunity interaction are focused on DDR deficiency including MMRd, homologous recombination - deficiency (HRD), and deleterious DDR mutations (23, 36). However, the role of overexpressed DDR profiles has been hardly investigated in immuno-oncology. We calculated the tumor-intrinsic signature scores of the urothelial cancer patients from a PD-L1 blockade cohort and the analysis indicated that DDR-related signatures were remarkably up-regulated in immunotherapy responders (Supplementary Figure 5A). Patients in response to anti-PD-L1 therapy exhibited higher signature scores of cell cycle, MMR, and homologous recombination signatures. Further, favorable therapy responses also correlated with increased DDR scores (Supplementary Figure 5B). Then the samples were divided into high and low DDR score groups according to the median score value. Intriguingly, high DDR score patients showed elevated TMB and TNB levels and prolonged survival (Supplementary Figures 5C–E). In summary, DDR score may predict the sensitivity of checkpoint immunotherapy in certain cancer types.