A DNA Repair and Cell-Cycle Gene Expression Signature in Primary and Recurrent Glioblastoma: Prognostic Value and Clinical Implications

Despite surgical resection and genotoxic treatment with ionizing radiation (IR) and the DNA alkylating agent temozolomide, glioblastoma (GBM) remains one of the most lethal cancers. Although occasional long-term survivors are reported (1), patients with GBM have a poor median survival (<1 year; refs. 2, 3), and all patients ultimately succumb due to treatment resistance and tumor relapse. Resistance to chemoradiation is promoted by complex DNA repair mechanisms, including O⁶-methylguanine-DNA methyltransferase (MGMT), which mediates the direct removal of O⁶-methylguanine (O⁶-meG), the most cytotoxic lesion induced by temozolomide. In the absence of MGMT, processing of O⁶-meG by the mismatch repair (MMR) pathway ultimately leads to perturbations of the replication fork and double-stranded DNA breaks (DSB) that require complex machineries for their repair. The other lesions induced by temozolomide are repaired mainly through base excision repair (BER) or direct removal mechanisms catalyzed by the DNA demethylases ALKBH2/3 (reviewed in ref. 4).

Concurrent and maintenance temozolomide was introduced to first-line treatment for GBMs in 2005 (2). Prior to this date, temozolomide and other alkylating agents were mainly used as second-line therapy in patients with recurrent tumors (5, 6). In the clinic, two populations of patients can be distinguished on the basis of methylation of the MGMT gene promoter-associated CpG island. Although patients with MGMT promoter–unmethylated GBM do not benefit from temozolomide, epigenetic silencing of MGMT (observed in about 40% of patients with GBM) confers a small but significant survival benefit (2.5 months) in patients exposed to temozolomide and IR compared with patients treated with IR only (7). However, improved therapeutic strategies are clearly needed in all cases.

Targeting components of the DNA damage response (DDR), including modulation of cell cycle and mitotic progression and genetic stability (8, 9), has emerged as an important therapeutic approach against many cancers. As a step toward improved strategies to undermine DNA repair or exploit specific vulnerabilities associated with GBM cells, we quantified the mRNA expression levels of a selection of genes covering the major DNA repair pathways, as well as important regulatory proteins and cell-cycle control genes, in GBM specimens and control, nontumor tissues. To this end, we exploited a cohort of paired GBM samples (i.e., matched primary and recurrent tumor from the same patients) that would allow us to address treatment-induced changes in gene expression, including biopsies obtained from patients treated with radiotherapy only (pre-2005) or with radiotherapy plus temozolomide (post-2005) so that the impact of temozolomide could be investigated.

The Köln cohort is described in Table 1. The expression data (Illumina microarray) pertaining to the 27 DNA repair and cell-cycle gene signature together with the relevant clinical data of the Heidelberg cohort (46 GBM pairs) are presented in Supplementary Table S1. Ethical guidelines were followed for patient sample collection and all samples were anonymized. Written informed consents were received from the patients, and the project was approved by local ethical committees in Cologne (Köln cohort, Application No. 03-170) or Heidelberg (Heidelberg cohort, Application No. 207/2005). Research was conducted according to the principles expressed in the Declaration of Helsinki. The Wang cohort is composed of 61 GBM pairs from three individual RNA-seq datasets (Korean_SMC, HF_MDA, TCGA_GBM) analyzed by RNA sequencing and described in ref. 10, and which are publicly available in the GlioVis database (http://recur.bioinfo.cnio.es).

All patients with GBM of the Köln cohort underwent surgical resection using standard craniotomy. The extent of resection was evaluated through pre- and postoperative MRI. Snap-frozen biopsies were used for analysis only if examination of their formalin-fixed paraffin-embedded counterpart by 2 independent neuro-histopathologists revealed >80% tumor cells. MGMT promoter methylation and IDH1 mutation status were assessed by DNA pyro-sequencing as described previously (11, 12). IHC stainings for Ki-67, PTTG1, AURKA, AURKB, and CENPA were performed on representative tissue sections from selected cases of primary and recurrent tumors from G1 and G3 groups. IHC was performed on an automated immunostainer (DAKO) using the UltraVision Quanto horseradish peroxidase detection system with 3,3′-diaminobenzidine tetrahydrochloride as chromogen (Thermo Fisher Scientific) and pretreatment of the sections by heating them in citrate buffer pH 6.0. The following primary mouse mAbs were used: anti-Ki67 (clone MIB-1, DAKO), anti-PTTG1 (clone DCS-280, MBL Co., 1:50), anti-AURKA (clone JLM28, Leica Biosystems, 1:50), anti-AURKB (clone mAbcam 3609, Abcam, 1:30), and anti-CENPA (clone 3-19, MBL Co., 1:100). The expression of Ki-67 was estimated using the following categories: <5%, 5%–10%, 10%–20%, 20%–30%, 30%–40%, 40%–50%, and >50% labeled nuclei.

Total RNA was extracted using the miRNeasy RNA Isolation Kit (Qiagen) according to the manufacturer's instructions. Samples with RNA integrity number (RIN) >7 and A_260/280 1.8–2.2 were sent to the VIB Nucleomics Core Facility for gene expression analysis using the direct multiplexed nCounter technology (NanoString; ref. 13) analysis and 154 gene probes (including 2 control, housekeeping genes) designed in silico. Gene expression data are available in the ArrayExpress database (https://www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-6425.

GBM patient–derived spheroid cell lines from G1 and G3 group tumors were obtained from the Heidelberg cohort, at the Department of Neurosurgery (University Clinic Heidelberg, Heidelberg, Germany). The cell lines were controlled for potential cross-contamination using short tandem repeat DNA typing and tested free of Mycoplasma. Cells were cultured in serum-free CSC medium (DMEM F-12 (Biowest L0093-500), supplemented with 20% BIT-100 (Provitro 204/3100), 1 U/mL Heparin (Sigma H3149-25KU), 4 mmol/L Ultraglutamine (Lonza BE17-605E/U1), 100U Pen-Strep (Lonza DE17-603E), 20 ng/mL EGF (Provitro 1325950500), and 20 ng/mL FGF (Miltenyi Biotec 130-093-841) in a humidified 37°C incubator under normal conditions (5% CO₂).

The G1 or G3 status of the two cell lines was ascertained by qRT-PCR analysis of the following genes from the signature: CCNA2, CDC25C, EME1, and TOP2A. To this end, RNA was prepared using TRIzol reagent (Ambion 15596018) followed by retrotranscription using the iScript cDNA Synthesis Kit protocol (Bio-rad 1708891). cDNA was then subjected to qPCR reaction in a ViiA7 Real-time PCR system (Applied Biosystems), using the FastSYBR green mix (Applied Biosystems 4385612) and the following primers: CCNA2 (Forward (F): AAGACGAGACGGGTTGC, Reverse (R): GGCTGTTTACTGTTTGCTTTCC), CDC25C (F: GACACCCAGAAGAGAATAATCATC, R: CGACACCTCAGCAACTCAG), EME1 (F: CTCCATGATACCCCAGAGAGG, R: CCTGGACCTTCTGACTCGG, TOP2A (F: ACAAGACATCAAAGTGAAGTAAAGCC, R: GCAGACTCAAAACACAGACAAAGC). The housekeeping genes used for normalization were GAPDH (F: CATGAGAAGTATGACAACAGCCT, R: AGTCCTTCCACGATACCAAAGT) and Ezrin (F: TGCCCCACGTCTGAGAATC, R: CGGCGCATATACAACTCATGG). Comparison of the relative expression of each gene in the G1 and G3 cell line was carried out using the 2^(−ΔΔC_t) method. Western blot analysis of Ki-67 expression levels in total cell extracts was carried out using the following mAbs: anti-Ki67 (clone MIB-1, DAKO, 1:100), anti-GAPDH (clone 5174S, Cell Signaling Technology, 1:1,000).

For cytotoxicity assays, 4,000 cells/well were seeded in triplicates in a 96-well plate format (Greiner 655101) followed by incubation with the selected drugs for 48 or 72 hours. Cytotoxicity was assayed using the WST-1 reagent (Roche 11644807001) and absorbances were determined using the Clariostar reader (BMG Labtech).

Differentially expressed genes were identified following raw RNA counts normalization, using the DESeq2 package (raw RNA counts; ref. 14) or limma package (15). Significant genes were identified using FDR and fold change. Functional annotation of the genes was carried out using the topGO R package. P values in the functional enrichement analysis were corrected for multiple hypothesis testing by FDR.

A thorough coexpression analysis was used to assess a gene expression pattern able to segregate the Köln cohort. First, we identified a cluster of genes with highly correlated profiles as determined by calculation of the coefficients of determination (R²) observed between log-transformed gene profiles. Second, this cluster of 52 highly correlated genes was considered for the classification of the clinical samples using nonnegative matrix factorization (NMF; ref. 16). After 1,000-fold NMF iterations with different initial estimations, the samples that clustered together in 90% of the runs were assigned to two independent groups (G1 or G3), whereas samples showing less concordance were assigned to an “undefined” G2 group. Refinement of the original 52-gene signature was achieved using the method developed by Heinäniemi and colleagues (17). In this analysis, gene pairs are randomly assigned, tested for specificity to our signature, and sorted according to their potential to participate in the patient clustering determination. This resulted in a refined signature with 27 genes that have the same capacity to segregate our cohort in 2 distinct groups as the initial 52-gene signature, as determined by NMF analysis. The robustness of the patient grouping obtained with both signatures was assessed by leave-one-out cross validation (LOOCV) and 5-fold cross validation. For the validation of the signature in the Heidelberg, Wang and Murat cohorts, and The Cancer Genome Atlas (TCGA) GBM RNA expression dataset, NMF was used to group GBM specimens based on the expression of the 27 gene signature, as described above.

Univariable Cox proportional hazards models were built using the R package survival (18), to investigate the potential impact of the following clinical factors on survival: gender, age, MGMT status and Karnofsky Performance Scale (KPS). Of these, age and KPS were considered as continuous variables, while others as categorical.

Correlation between the gene expression levels pertaining to the 27 gene signature among the Köln cohort and the validation datasets (Heidelberg and Wang cohorts) was assessed by calculating the logFC for each gene, followed by linear regression analysis.

To investigate DNA repair and cell-cycle control gene expression in GBM, we used the nCounter technology, a robust and sensitive method allowing multiplexed analysis of a panel of selected genes (19), to analyze RNA extracted from the Köln cohort (Table 1) composed of samples from paired GBM biopsies (n = 66, of which, 18 were from temozolomide-naïve patients and 48 from temozolomide-treated patients) as well as unpaired biopsies (n = 12) and control, tumor-adjacent tissues (n = 9). All patients received radiotherapy and harbored wild-type IDH1, as determined by DNA sequencing. Because of the nCounter gene expression format, we focused on 154 genes encompassing the major DNA repair pathways: base excision repair (BER), nucleotide excision repair (NER), mismatch repair (MMR), homologous recombination (HR), nonhomologous end joining (NHEJ), direct repair enzymes (e.g., MGMT), Fanconi anemia (FA), as well as a selection of genes encoding effectors and regulators of cell cycle, DNA replication, DDR, centromere, and centrosome dynamics (Supplementary Table S2). These genes were chosen based on the available literature on the response of GBM cells to chemoradiation and in particular to temozolomide (4).

In preliminary comparisons of tumor samples with control, tumor-adjacent tissues, detailed in Supplementary Data S1, we identified lists of differentially expressed genes (DEG) associated with primary (72 genes) and recurrent GBM tumors (71 genes). These lists displayed extensive overlap, with 67 genes in common that exhibited comparable fold change direction and intensities. Accordingly, no DEGs were identified upon direct comparison of primary and recurrent GBM samples. Notably, no temozolomide-associated changes in gene expression were observed when temozolomide-naïve and temozolomide-treated recurrent biopsies were compared.

Failure to identify temozolomide-associated gene expression changes using global DEG analysis of unstratified specimens prompted us to examine whether stratification would help uncover such genes in a subset of patients. To address this issue and with the aim of relating our findings to tumor recurrence, we first subjected the gene expression data from the paired biopsies (66 specimens) to a coexpression analysis and obtained a set of 52 genes associated with highly correlated profiles (Supplementary Fig. S1A). We then used this set of genes for the stratification of our specimens using NMF. The resulting heatmap (Supplementary Fig. S1B) revealed that this 52-gene signature segregated the biopsies into 2 well-defined groups (hereafter called G1 and G3) that displayed an inverse expression pattern of the gene signature, leaving 10 samples with a more neutral profile in a separate group called G2. Our 52-gene signature was robust, as evaluated by 2 forms of cross-validation [mean sample misclassification error: 0.0606 (LOOVC), 0.087 (5-fold cross-validation)].

The signature contained 2 subsets of genes [component A (n = 27) and B (n = 25) of Supplementary Fig. S1C] characterized by inverse expression in G1 and G3. Furthermore, component A was found to contain all the cell cycle–associated genes in the signature.

Inspired by Heinäniemi and colleagues (17), we refined the 52-gene signature by analyzing the predictive power of all gene pairs able to reproduce the initial classification. We identified a subset of 27 genes within which each gene pair displayed high prediction ability [mean misclassification error: 0.0758 (LOOCV), 0.086 (5-fold cross-validation)] and, interestingly, highly positive coexpression (Fig. 1A). Notably, 21 of those 27 genes belonged to component A (Supplementary Fig. S1C) of the original signature. The resulting heatmap generated by NMF clustering is presented in Fig. 1B, again revealing G1 and G3 groups presenting an inverse gene expression pattern, and a third less-defined group (G2), as observed with the original signature (Supplementary Fig. S1B).

All 3 groups contained both primary and recurrent samples. Indeed, stratification into 3 groups was also obtained when primary and recurrent biopsies were considered separately, and when an independent cohort of primary GBMs (20) was examined (see below). However, whereas some patients (13/33) had both their primary and recurrent biopsies in the same group, in 20 of 33 patients, the recurrence was traced to a group distinct from that of the primary tumor (Fig. 1C), a notion hereafter referred to as group migration. These observations indicate that relapse was often associated with significant alterations in DNA repair and cell-cycle gene expression. The degree of alteration culminated in patients whose paired biopsies showed inverse gene expression of the signature (Fig. 1B). Similar group migrations were observed with the original 52-gene signature (Supplementary Fig. S1D).

Our 27-gene signature contains important effectors and regulators of the cell cycle, centromere and centrosome dynamics, chromosome segregation, and mitosis (AURKA, AURKB, CCNA2, CCNB1, CDC25C, CDC6, CDK1, CENPA, CENPF, MKI67, PCLAF, PLK1, PTTG1, TOP2A), as well as genes encoding crucial HR factors such as the RAD51 recombinase, the chromatin remodelers RAD54B and RAD54L, and enzymes involved in Holliday junction resolution (EME1/MUS81 complex) and/or NER (ERCC3(XPB), ERCC4(XPF)). Also in the signature were genes encoding the DNA glycosylase NEIL3, Fanconi Anemia factors (FANCD2, UBE2T), the ubiquitin protein ligase UBE3B and 2 specialized DNA polymerases, POLM, and POLQ, involved in NHEJ pathways of DSB repair. None of these genes have been associated to sites of frequent copy number alterations (21) or found to be significantly mutated (22) in GBM.

Given the inverse gene expression profile displayed by the G1 and G3 groups and the large number of genes involved in cell cycle and mitosis regulation in our signature, including MKI67 encoding the proliferation marker Ki-67, we further investigated the expression of this marker at the protein level. IHC staining of Ki-67 revealed that 95% of the G1 samples exhibited a low Ki-67 index (<10%), whereas 72% of G3 samples displayed a high Ki-67 index (>30%; Fig. 2), in line with the gene expression data and suggesting that the net effect of the observed upregulation of the cell-cycle genes in G3 was increased proliferation. Additional immunostainings for PTTG1, AURKA, AURKB, and CENP-A were carried out using selected primary and recurrent tumors from the G1 and G3 groups (Fig. 2), globally validating the expression of these genes at the protein level.

We next sought to validate our signature by challenging two datasets: the Heidelberg cohort (biopsy pairs n = 46) and a heterogenous collection described in ref. 10, hereby referred to as the Wang cohort (biopsy pairs n = 61).

NMF clustering of the Heidelberg cohort generated 3 groups that exhibited profiles reminiscent of the G1, G2, and G3 groups (Supplementary Fig. S2A) and displayed significant group migration (Supplementary Fig. S2B). Similar observations were made with the Wang dataset (Supplementary Fig. S2C and S2D). Furthermore, although the heatmaps generated by NMF suggested that, when considered individually, not all the genes of our signature performed equally well in the validation datasets, the behavior of the gene signature was concordant among all datasets, as revealed by linear regression analysis (P values: 0.0019 for Köln and Heidelberg; 6.1 × 10⁻⁶, for Köln and Wang; Supplementary Fig. S2E and S2F). Thus, our 27-gene signature allowed classification of GBM tumors based on DNA repair and cell-cycle gene expression, and performed equally well with datasets generated by different RNA analysis platforms.

We next asked whether the stratification afforded by our signature now allowed temozolomide-associated changes to be uncovered in specific subsets of samples. Comparable numbers of recurrences from temozolomide-naïve and temozolomide-treated patients were found within each group (Fig. 1B), suggesting that the clustering was independent of temozolomide treatment. Moreover, within-group comparisons of the recurrences from temozolomide-naïve and temozolomide-treated patients did not uncover DEGs among the whole gene set, indicating that the initial failure to observe temozolomide-associated genes in our cohort was not due to lack of patient stratification. Finally, when we considered the 31 of 33 patients with known MGMT methylation status of the primary tumor, we found no statistically significant differences in the representation of MGMT-methylated and -unmethylated specimens among the 3 groups. The lack of annotations on MGMT promoter methylation status in the Heidelberg and Wang cohorts led us to ascertain this observation using the cohort of primary GBM specimens with well-defined MGMT status described by Murat and colleagues (20). NMF clustering of the 84 samples of this cohort using our 27-gene signature again identified 3 groups, including well-defined G1 and G3 groups (Supplementary Fig. S3A). In these 2 groups, the relative distribution of specimens with methylated and unmethylated MGMT promoter was identical (Supplementary Fig. S3A), indicating that the clustering was independent of MGMT promoter methylation status.

We next examined progression-free survival (PFS) and overall survival (OS) in the 3 groups. When only primary tumors were considered, analysis of these endpoints showed no significant differences between the groups, neither in the paired GBM cohorts, nor in two independent cohorts of primary GBMs [Murat and colleagues (Supplementary Fig. S3B; ref. 20)], and the TCGA GBM RNA expression dataset (Supplementary Fig. S4). In contrast, when only the recurrent biopsies of the Köln cohort were analyzed, we found a significant difference in PFS between patients with G1 and G3 biopsies: PFS was longer among patients with a G3 recurrence (12.3 months), compared with those with a G1 recurrence (6.05 months, Fig. 3A). A similar trend was observed with OS, although not statistically significant (inset in Fig. 3A). Importantly, univariable Cox analysis of the Köln dataset interrogating age, gender, MGMT promoter status, KPS, and grouping identified the recurrence association with a specific group as the sole parameter influencing PFS (P = 0.019; Supplementary Table S3).

Group migration analysis in the Köln dataset revealed that all transitions were possible (Fig. 1C), and this was corroborated in the other datasets (Supplementary Fig. S2B and S2E). To further explore the impact of group migration on the survival parameters, we assembled patients of the Köln cohort based on the migration behavior of their recurrence, thus creating 4 categories (A–D) corresponding, respectively, to those patients with altered recurrence in G1 (A), G2 (B), G3 (C) or with unaltered recurrence (D). Pairwise analyses revealed that the worst PFS and OS were associated with category A (i.e, progression to G1; Fig. 3B). These findings were reinforced by the observation that both PFS (Fig. 3C) and OS (Fig. 3D) were significantly worse in patients with G1 recurrences (rG1) originating from G3 primaries (pG3; defined as pG3-rG1, red line) than in patients with G3 recurrences originating from either G1 or G3 primaries (pG1-rG3, black line and pG3-rG3, green line), suggesting that these patients contributed the major determinant behind the poor prognosis associated with category A. Taken together, our data suggest that alterations of the DNA repair and cell-cycle gene expression signature at relapse are associated with significant changes in prognosis, with a transition to G1 translating into the poorest prognosis. Finally, the G1 and G3 groups presented similar proportions of recurrent tumors from temozolomide-treated and -naïve patients (Fig. 1B), indicating that the observed difference in survival between G1 and G3 was not related to temozolomide treatment.

That a G1-type expression pattern at relapse was associated with worse survival parameters was also observed in the Heidelberg cohort where statistically significant differences in PFS and OS were seen both when all recurrences were considered (Fig. 3E) and when group migration was highlighted (Fig. 3F).

Having shown that a recurrent G1 tumor was associated with worse prognosis (Fig. 3C, D, and F), we next examined whether it was the migration toward the G1 group or the mere presence in G1 that was associated with poor prognosis, by comparing pG1-rG1 and pG2/G3-rG1 patients. As only one patient of the Köln cohort belonged to the first category, we focused on the Heidelberg cohort. Although pG1-rG1 patients (n = 3) had worse prognosis than pG2/G3-rG1 patients (n = 14; PFS: 6 vs. 15.5 months; OS: 17 vs. 25.4 months), the differences were not statistically significant.

No association could be established between gene expression at relapse and survival parameters in the Wang cohort. Based on published work (23), we considered differences in the age distribution as a potential explanation for this finding. Indeed, while there was no statistical difference in age between the Köln cohort (median age: 58 years, range: 30–73) and the Heidelberg cohort (59 years, range: 32–81; P = 0.73, Wilcoxon test), the difference was significant between the Köln and Wang cohorts (median age: 52 years, range: 29–74; P = 1.2 × 10⁻⁵, Wilcoxon test), with the age distribution in the Wang cohort being shifted towards younger patients. In addition, unlike the Köln and Heidelberg cohorts, Cox regression analysis indicated that age had an impact on survival of patients in the Wang cohort (P = 1.04 × 10⁻⁶). Age is a strong prognostic factor among patients with GBM, with younger age correlating with improved survival (24). Thus, although the contribution of other clinical factors cannot be excluded, our observations suggest that our failure to validate the prognostic impact of our signature in the Wang cohort may in part be due to the increased proportion of younger patients in this cohort, which obscures the impact of our patient stratification.

Recent refinement of the gene expression subtypes associated with GBMs has led to three subtypes: proneural, classical, and mesenchymal (10). Having validated our signature using cohorts from this study, we used its annotations to probe the association between our groups and the transcriptional subtypes. Notably, when all the biopsies were considered (i.e., primary and recurrent specimens), we observed specific enrichment of mesenchymal-, classical- and proneural-type tumors in group G1 (49%, P = 0.06), G2 (52%, P = 0.01), and G3 (42%, P = 0.01), respectively (Supplementary Fig. S2D).

We next investigated the biological processes that distinguished the G1 and G3 groups displaying inverse expression profiles of our gene signature, using the RNA expression dataset from the Wang cohort whose stratification is shown in Supplementary Fig. S2. We identified 2061 DEGs between the two groups (FDR < 0.001; |logFC| > 1), of which, 253 were upregulated in G3 and 1808 in G1 (Supplementary Fig. S5A; Supplementary Table S4). Gene Ontology (GO) enrichment analysis (Supplementary Fig. S5B) indicated that GO terms corresponding to biological processes related to cell cycle and cell division, chromosome organization and segregation, DNA replication and repair, and chromatin assembly and dynamics predominated among the DEGs upregulated in the G3 group. These results are in full agreement with the initial distinction of the G1 and G3 groups obtained with our gene signature. In contrast, the list of biological processes significantly enriched among the DEGs upregulated in the G1 group included terms related to vesicle-mediated transport, apoptosis and autophagy, cellular adhesion, response to (chemical) stimulus, and response to stress. These results are consistent with the notion that G3 group tumors develop cellular programmes to sustain high proliferation, in contrast to G1 tumors, which engage critical survival cellular processes.

As a step to determine whether specific genetic aberrations were associated with tumor development and/or progression within the groups identified by our signature, we took advantage of the availability of exome sequencing and copy number data for a subset of paired tumors of the Wang cohort characterized by Kim and colleagues (25). We generated an oncoprint where the grouping assigned by NMF based on our signature was superimposed on the panel of genetic aberrations related to frequently amplified, deleted and mutated gene components of the major signaling pathways involved in the pathogenesis of GBM, previously determined for these samples (Fig. 4A). Notable differences between the G1 and G3 group tumors are schematically summarized in Fig. 4B. Unlike G3 tumors, all G1 tumors displayed deletion of PTEN and also a high frequency of EGFR and PDGFRA amplification, strongly implicating growth factor tyrosine kinase receptor pathways and the PI3K/PTEN/AKT pathway in the etiology of these tumors (26). In contrast, G3 tumors were characterized by an increased number of cases with RB1 deletion, as well as a high frequency of mutations in TP53, PTEN, and NF1, suggesting that PI3K/PTEN/AKT, RB/CDKN2A-p16, and TP53 pathways contributed to the development of G3 tumors. Tumors from the G2 groups presented a composite pattern, for instance displaying increased frequency of TP53 mutations compared with G1 tumors, but also higher frequency of PTEN deletion compared to G3 tumors (Fig. 4A). Finally, when group migrations were considered in this small cohort, the majority of cases did not associate with mutations/alterations in the genes characterized by Kim and colleagues (25). In addition, no group migration involving pG3-rG1 was observed in this cohort. However, we noted that one case of pG1-rG3 migration involved transitions in the clonal status of TP53 and PTEN mutations, from subclonal in the initial G1 tumor to clonal in the recurrent G3 tumor.

We next wanted to know whether our 52- and 27-gene signatures exposed group-specific vulnerabilities against clinically relevant genotoxic agents as well as inhibitors of the DDR and perturbators of the cell cycle. Indeed, increased expression of TOP2A in the G3 group suggested that G3 tumors might be more sensitive to topoisomerase II inhibitors such as etoposide (27). Likewise, the downregulation of NER genes ERCC3/XPB and ERCC4/XPF in this group suggested that agents such as cisplatin that cause interstrand crosslinks requiring NER for their repair (28, 29), might be more efficient against G3 tumors. We also considered inhibition of RAD51, involved in the recombinational repair of DSBs (30), which was upregulated in group G3 compared to G1, as well as inhibitors of mitotic kinases, given the strong expression of AURKA, AURKB, CDK1, and PLK1 in G3.

Having successfully established two cell lines, NCH741f and NCH481, from clinical biopsies obtained from G1- and G3-group patients, respectively, of the Heidelberg cohort, we first verified their G1 and G3 status by RT-qPCR analysis of selected genes from the 27-gene signature, including CCNA2 and CDC25C, which together form a gene pair with high predictive ability (mean misclassification error: LOOCV = 0.07576; 5-fold cross validation = 0.09774 ± 0.0015; Fig. 5A, left). We also performed Western blot analysis of Ki-67 showing increased Ki-67 levels in the G3 cell line compared with the G1 cell line (Fig. 5A, middle), consistent with the calculated doubling times for these cell lines (NCH481: ca. 57 hours; NCH741f: ca. 104 hours). Importantly, these cell lines were obtained and maintained in serum-free medium conditions, which is crucial to preserve the biological status and behavior of GBM cells, including their response to therapeutics (31). Under these conditions, the cells grew as spheroids (Fig. 5A, right).

We next treated the 2 cell lines with the DNA-damaging agents etoposide and cisplatin, as well as inhibitors of RAD51 (i.e., RI-1) and mitotic kinases [i.e., anti-AURKA (alisertib) and pan-AURK (tozasertib)] to which we expected cells from the G3 group to be more sensitive compared with G1 cells. Given that inhibitors of the DNA repair factor PARP (PARPi; ref. 9) are currently in clinical trials against GBM (32), we also tested the PARPi olaparib. In this case, however, we could not predict differences in sensitivity to the PARPi among G1- and G3-group tumor cells based on the expression of our gene signature components because synthetic lethality with PARP inhibition can be mediated by various defects in a plethora of DNA repair pathways (33, 34). The data from WST-1 cytotoxicity assays illustrated in Fig. 5B, show the statistically significant increase in sensitivity of G3-group derived cells relative to G1-group cells, to all tested compounds with the exception of olaparib, for which no significant differences were observed, in line with our gene expression data.

Taken together, these data provide experimental evidence to suggest that analysis of GBM biopsies using our DNA repair and cell-cycle gene signature may help identify novel therapeutic strategies against subsets of GBM tumors.

Previous studies sought to identify molecular profiles specific to resistance to chemoradiotherapy by interrogating clinical outcome data and gene expression data generated from primary GBM specimens (i.e., collected before treatment; ref. 20). In this study, we focused on a cohort of paired GBM samples from patients treated with radiotherapy or radiotherapy plus temozolomide with the aim of exploring alterations in DNA repair and cell-cycle gene expression in primary GBMs and their recurrences. Although the use of paired biopsies is gaining momentum, their availability is still limited by the speed of progress in the biobanking of reliable and clinically annotated specimens. For GBM, paired biopsies predating the introduction of temozolomide in the clinic represent a scarce resource.

Our analysis did not reveal temozolomide-related alterations in the studied genes, or expression patterns predictive of response to temozolomide. This observation suggests that the dynamics of mRNA expression in response to temozolomide treatment cannot be assessed in clinical samples collected at the time of second surgery, most likely due to the delay between the cycles of temozolomide administration, tumor relapse, and sampling of the recurrent biopsy. Experiments to investigate this question would benefit from patient-derived xenograft models where time courses of gene expression in response to treatment can be carried out, as underlined by studies in ovarian (35) or breast (36) cancers. The mutagenic impact of temozolomide as a driving mechanism in tumor progression and the selection of temozolomide-resistant clones through defective MMR is well documented in gliomas (37–39) and GBM (40, 41). Future work should focus on the mutational status of the MMR genes and other mutations induced by temozolomide in our cohorts.

Our 27-gene signature clustered primary and recurrent GBMs according to DNA repair and cell-cycle gene expression. Performance comparison with other signatures is hampered by the current paucity of studies on paired biospies. Thus, a 412-gene classifier has been described that segregated paired GBM specimens into 2 groups, revealing cases of group migration; however, no association between patient grouping and survival was reported (42). Although our gene signature could not be used to predict relapse-free survival based on the profiling of primary GBM samples, its expression pattern at the time of recurrence correlated with prognosis. As our signature contains 14 cell cycle genes that display high positive coexpression, it is striking that decreased expression of these genes at relapse in G1 was associated with shorter PFS. It is possible that the high Ki-67 proliferation index seen in the G3 recurrent tumors is associated with exacerbated genetic and/or chromosomal instability, leading to increased mitotic catastrophe and cell death. Increased chromosomal instability has been described in recurrent GBM (25, 43). In addition, deregulation of CENPA and CENPF, present in our signature, could contribute to the observed phenotype, because misregulation of centromere and kinetochore function leading to chromosomal instability contributes to tumor development and response to chemoradiation (44). Alterations in the division mode and differentiation status of the tumor cells and their subpopulations may also affect chromosomal instability and proliferation in ways that remain to be explored. As no standard second-line treatment of GBM has yet been determined, patients with recurrent GBM often receive complex modalities (5, 6). Because the survival outcome of these patients may benefit from second-line therapy, in depth analysis of these modalities may help determine the extent to which the improved OS of patients with a G3-group recurrence reflects the impact of DNA repair and cell-cycle gene expression on the tumor cell response to second-line treatment. Such an analysis was not possible in our cohorts, due to the lack of documented patient records and the large panel of drugs that were considered as second-line treatment, resulting in multiple subgroups of patients whose small size precluded any statistically significant analysis.

It is notable that the G1 group from the Wang cohort was enriched in mesenchymal-type specimens, since this subtype has been shown to be associated with worse PFS in recurrent GBM (ref. 45 and references therein). Furthermore, analysis of key molecular alterations that underlie GBM in a small subset of tumors from this cohort suggests that there are distinct molecular drivers behind the groups identified by our signature. However, extensive molecular characterization of the paired biopsies of our cohorts, for example, by large-scale next-generation sequencing approaches, will be required to understand the full impact of chemoradiation on genomic alterations and reveal the molecular drivers behind the groups detected in our studies and their associated prognosis. Such comprehensive analyses are currently pursued in large collaborative efforts such as the one undertaken by the Glioma Longitudinal Analysis (GLASS) consortium (46).

Our gene signature exposes vulnerabilities against specific DDR inhibitors, perturbators of mitosis, and/or genotoxic agents, including several mainstream chemotherapeutics (etoposide, cisplatin, carboplatin) previously considered, either alone (47–49) or in combination (50, 51) against GBM. It is tempting to speculate that the heterogeneity of the cohorts used in these studies contributed to the limited efficiency attributed to these compounds and that future studies might benefit from patient's tumor stratification based on our signature. Potential modalities also involve the targeting of crucial components of the DNA damage response including RAD51 (9), which our data suggest could prove an attractive target to undermine HR-mediated repair of replication fork-associated DSBs in highly proliferative cells (e.g., in G3-group tumors) or to sensitize these cells to IR and genotoxic agents known to induce DSBs. The identification of patients that are likely to benefit from PARP inhibition remains a challenge (33, 34), and we could not predict differences in sensitivity to the PARPi olaparib used as a single agent based on the expression of our gene signature components. Future experiments will need to consider PARPi in combination with DNA-damaging agents or other DDRi. Our data also lend support to the use of mitosis perturbators against GBM. Indeed, our signature contains crucial regulators of this process, including PTTG1/securin and the mitotic kinases AURKA, AURKB, CDK1 and PLK1, which have gained recent attention as therapeutical targets in GBM. Thus, in vitro and xenograft animal model studies have advocated the use of AURK inhibitors in combination with radiotherapy and temozolomide against GBM cells (52–55). Our data suggest that, when used as single agents, mitotic kinase inhibitors may be particularly effective against G3-group tumors. Finally, although our experimental data, obtained with genotoxicants and inhibitors used in monotherapy formats, suggest that G3-group tumor cells are particularly amenable to DNA repair–based strategies, future experiments using combination therapy may unveil modalities targeting specifically G1-group cells. How our gene signature may contribute to precision medicine for GBM is summarized in Supplementary Fig. S6. We are aware that our stratification strategy still leaves aside a considerable subset of patients (e.g., from the G2 group). However, in view of the fact that GBM remains an incurable disease, our results may carry significant clinical prospects in the fight against GBM. The modalities discussed here in the framework of our gene signature might be applicable for the management of both primary and recurring GBM.

G. Reifenberger reports receiving commercial research grants from Merck and Roche, has received speakers bureau honoraria from Abbvie, and is a consultant/advisory board member for Abbvie. A.J. Chalmers reports receiving commercial research grants from AstraZeneca and Vertex Pharmaceuticals, has received speakers bureau honoraria from Tesaro and Bayer, and is a consultant/advisory board member for Hox Therapeutics and AstraZeneca. No potential conflicts of interest were disclosed by the other authors.

Conception and design: A.J. Chalmers, R. Goldbrunner, E. Van Dyck

Development of methodology: M. Gobin, P.V. Nazarov, E. Van Dyck

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Gobin, R. Warta, M. Timmer, G. Reifenberger, J. Felsberg, A.J. Chalmers, C.C. Herold-Mende, R. Goldbrunner, S.P. Niclou, E. Van Dyck

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Gobin, P.V. Nazarov, R. Warta, G. Reifenberger, J. Felsberg, A.J. Chalmers, C.C. Herold-Mende, E. Van Dyck

Writing, review, and/or revision of the manuscript: M. Gobin, P.V. Nazarov, R. Warta, M. Timmer, G. Reifenberger, J. Felsberg, L. Vallar, A.J. Chalmers, C.C. Herold-Mende, E. Van Dyck

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Timmer, S.P. Niclou, E. Van Dyck

Study supervision: S.P. Niclou, E. Van Dyck

This work was supported by Télévie-FNRS (grant 7.6532.16 to M. Gobin and E. Van Dyck), the Luxembourg National Research Fund (FNR grant C17/BM/11664971/DEMICS to P.V. Nazarov), and the Anni Hofmann Stiftung (R. Warta).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.