Association of PARP1 polymorphisms with response to chemotherapy in patients with high‐risk neuroblastoma

Abstract The genetic aetiology and the molecular mechanisms that characterize high‐risk neuroblastoma are still little understood. The majority of high‐risk neuroblastoma patients do not take advantage of current induction therapy. So far, one of the main reasons liable for cancer therapeutic failure is the acquisition of resistance to cytotoxic anticancer drugs, because of the DNA repair system of tumour cells. PARP1 is one of the main DNA damage sensors involved in the DNA repair system and genomic stability. We observed that high PARP1 mRNA level is associated with unfavourable prognosis in 3 public gene expression NB patients’ datasets and in 20 neuroblastomas analysed by qRT‐PCR. Among 4983 SNPs in PARP1, we selected two potential functional SNPs. We investigated the association of rs907187, in PARP1 promoter, and rs2048426 in non‐coding region with response chemotherapy in 121 Italian patients with high‐risk NB. Results showed that minor G allele of rs907187 associated with induction response of patients (P = .02) and with decrease PARP1 mRNA levels in NB cell line (P = .003). Furthermore, rs907187 was predicted to alter the binding site of E2F1 transcription factor. Specifically, allele G had low binding affinity with E2F1 whose expression positively correlates with PARP1 expression and associated with poor prognosis of patients with NB. By contrast, we did not find genetic association for the SNP rs2048426. These data reveal rs907187 as a novel potential risk variant associated with the failure of induction therapy for high‐risk NB.


| INTRODUC TI ON
Neuroblastoma (NB) is the most frequent malignant tumour in paediatric age arising from neural crest cells, precursors of the sympathetic nervous system. 1 It is an enigmatic tumour due to its high genetic heterogeneity and because of the complexity to develop a successful therapy by clinicians and researchers. The therapeutic regimens for high-risk patients (defined following the International Neuroblastoma Risk Group (INRG) classification system 2 ) include in the first phase an induction therapy with an intensive cycle of chemotherapeutic agents. Afterwards, patients follow a consolidation phase with myeloablative therapy, stem cell transplantation (SCT) and radiation therapy to primary tumour and residual metastatic sites, followed by maintenance aimed at controlling minimal residual disease. 3,4 Analyses of recent published high-risk NB trials emphasize the large proportion of patients who do not continue beyond induction chemotherapy owing to inadequate response. For example, in the recent Children's Oncology Group (COG) trial examining purged versus non-purged peripheral blood stem cell transplantation, 25% of the 495 registered patients did not continue to treatment beyond induction chemotherapy, most commonly because of progressive disease (56%). 5 A larger proportion of patients (51.5% of 1,231) failed to continue beyond induction therapy in the International Society of Pediatric Oncology Europe Neuroblastoma (SIOPEN) high-risk NBL-1 trial. 6 Though the rates of patients continuing to consolidation therapy differ, partly owing to differing consolidation criteria the underlying issue remains that a substantial number of patients fail to respond to current high-risk NB induction therapy. Therefore, advanced NB remains one of the most unmanageable paediatric cancers with long-term survival rate below 50%. 7 The development of reliable biomarkers for clinical implementation is therefore a priority. MYCN amplification and segmental chromosomal aberrations are, so far, the most reliable genomic biomarkers for the patients' stratification and outcome prediction.
Recently, genome-wide association studies and high-throughput sequencing-based studies have highlighted that multiple DNA polymorphisms influence NB susceptibility and clinical phenotype [8][9][10][11][12][13] and that recurrent mutations of single genes are infrequent in primary NB with activating mutations in ALK and inactivating mutations in ATRX, and TERT rearrangements being the most frequent. [14][15][16][17][18] Gene expression-based studies suggest that among high-risk patients, gene signatures can identify children with higher risk disease who would benefit from new and more aggressive therapeutic approaches. [19][20][21] Despite these large efforts made to find genomic biomarkers for improving high-risk patient outcome, so far no study has searched for heritable variations able to predict the primary effect of chemotherapy.
One of the main reasons responsible for cancer therapeutic failure is the acquisition of resistance phenotypes to cytotoxic anticancer drugs. This is mainly because of the efficiency of the DNA repair system of cancer cells, which enhances the tolerance to DNA damages induced by chemotherapy and radiotherapy. 22,23 DNA damaging cancer therapeutics take advantage of overlapping DNA repair pathways, including base excision repair (BER), nucleotide excision repair (NER), double-strand break repair (DSBR) and mismatch repair (MMR) pathways. 24 As BER is one of the major DNA repair pathways, reducing BER capacity is a useful approach for cancer treatment. 25 PARP1 belongs to the family of the poly (adenosine diphosphate-ribose) polymerase (PARP) proteins, which are DNA damage sensors, with the ability to signal to downstream effectors and with that directly involved in genomic stability, DNA repair and apoptosis. 24 The roles of PARP1 in the DNA damage response have been studied extensively. Induction of various kinds of DNA damage results in rapid recruitment of PARP1 to sites of damage through its DNA-binding ability. 26 It is involved in (DSBs) in both homology-directed repair (HDR) and non-homologous end-joining (NHEJ) pathways. In addition, single-strand breaks (SSBs) are very rapidly detected and bound by PARP1, which also represents one of the components of the BER process simplifying the subsequent recruitment of BER proteins. 27 PARP inhibitors became interesting tools to boost the activity of cancer chemotherapy. Therefore, many clinical trials are examining their efficiency in combined therapy approaches in different sets of patients with cancer. 28 Taking into account of the evident number of patients with NB that do not reply therapies, there is interest to early identify the non-responder patients to allow their enrolment to suitable treatments. The purpose of our study is to identify functional single nucleotide polymorphisms (SNPs) of PARP1 able to predict the response to current induction therapy in patients with high-risk NB.

| Cataloguing of functional SNPs in PARP1
To identify functional SNPs, we listed 4983 SNPs in PARP1, within 1 ± Mb surrounding the gene and with minor allele frequency (MAF) greater than 10%. Thus, in order to identify SNPs in PARP1 gene that may be associated with NB patients' induction response, we performed a filtering strategy of PARP1 variants ( Figure 2 and Table S1).

| SNP genotyping
The SNPs rs907187 and rs2048426 were genotyped by TaqMan SNP Genotyping Assay as previously described 29 on 7900HT Realtime PCR System (Applied Biosystems). To monitor quality control, three DNA samples per genotype were genotyped by Sanger sequencing (3730 DNA analyzer; Applied Biosystems) and included in each 384-well reaction plate; genotype concordance was 100%.
To confirm genotypes, we sequenced 20 samples chosen randomly from responders and non-responders; concordance between genotype was 100%. Primer sequences are available upon demand.

| Construction of luciferase reporter gene plasmids
The genomic region of 1111 bp expanding from 555 bp upstream to 555 bp downstream the variant rs907187 was cloned upstream of the firefly luciferase gene. PCR primers contained recognition sites for NheI in the forward and XhoI in the reverse primer were designed to amplify 1153 bp from the genomic DNA of cell lines homozygous for the rs907187-G allele. After cutting the fragment with NheI and XhoI restriction enzyme (Biolabs), we cloned it into the pGL3-Basic Vector (Promega). The resulting plasmid containing the rs907187-G allele was site-specifically mutated to C alleles using Site-Directed Mutagenesis Kit (Stratagene). The sequence of each construct was confirmed by direct sequencing.

| In vitro functional analysis
SKNBE2 cells were transfected X-tremeGENE (Roche) with 1 ug of pGL3-Basic Vector rs907187-C and rs907187-G constructs. Cells were subsequently starved in serum-free medium for 8 hours.    The experiments were carried out in triplicate for each data point.

| Western blotting
The housekeeping gene β-actin was used as internal control.
Relative gene expression was calculated using the 2 −ΔCT method, where the ΔCT was calculated using the differences in the mean CT between the selected genes and the internal control (β-actin).

| Association between clinical outcome and PARP1 and PARP2 expression in NB patients
To make evident the association of PARP1 with clinical outcome of patients with NB, we evaluated its expression in three public gene expression datasets of NB (described in materials and methods). Data in Figure 1A; Figures

| Identification of functional SNPs in PARP1
To define a set of credible risk variants in PARP1 gene, that may be associated with NB patients' induction response, we selected 176 variants (Figure 2 and Table S1) Table S1) that altered the binding sites (prediction made by motif braker 31 ) of transcription factors experimentally (ENCODE project) found to occupy the same sites ( Figure 2; Table   S1). As last step, among the selected 9 variants, we considered those which affect the quantitative trait loci expression (eQTLs) of PARP1 (GTEx portal data in Table S2). This resulted in 3 highly significant polymorphisms in non-coding regions: rs892348, rs2048426 and rs907187 (coloured in yellow in the Table S1). We also selected 1 coding SNP (rs1136410) predicted highly pathogenic through CADD score (=32) (Figure 2 and Table S3). As reported in Figure S4, rs892348, rs907187 and rs1136410 are in strong LD (>0.9), so we decided to analyse rs907187, in PARP1 promoter, and rs2048426 in non-coding region.

| Association of rs907187 G allele with response to induction therapy
In the analysis of all 121 patients, the minor G allele of rs907187 (in the promoter region of PARP1) associated with a better response (P = .02, Table 2), by contrast, we found no genetic association for the SNP rs2048426 ( Table 2).
The multivariate analysis including the prognostic factors MYCN, age at diagnosis and 1p36 deletion confirmed the association between G allele and good response to the therapy (Table S4). We also tested whether the SNP genotypes were associated with overall and event-free survival, and prognostic factors ( Figure S5 and Table S5).
No significant associations were found.
We corroborated these data analysing the gene expression variation, using genome-wide expression and SNP arrays of NB tumours, demonstrating that the SNP affects the expression of PARP1.
In particular, the presence of the protective allele G correlated with decreased PARP1 mRNA expression in a set of 17 NB cell lines ( Figure 3A; P = .0003). This result was furthermore explored through SNP-expression correlation on Genotype-Tissue Expression (GTEx) project available through the GTEx portal. Notable, G allele of rs907187 correlated with a lower PARP1 expression in the analysis of 256 Nerve Tibia Tissues ( Figure 3B; P = .000012). Finally, to further assess the impacts of rs907187 on PARP1 expression, we performed a luciferase report gene assay and observed that the induction of promoter activity of the construct containing rs907187-G alleles was lower than that of the construct containing C alleles in NB SKNBE2 cells ( Figure 3C).
Together, these findings indicated that the decrease in PARP1 expression because of the rs907187-G allele may predispose NB patients to better response to induction therapy and support a potential role of PARP1 as a candidate gene in therapeutic failure of NB treatment.

| rs907187 is predicted to alter the binding site of E2F1
The tool TRAP 33 predicted a high binding affinity of the risk allele C to E2F1 (E2F transcription factor 1) that was not observed with G allele ( Figure 4A). The E2F1 transcription factor plays a crucial role in the control of cell cycle acting as an oncogene, 34-36 so we tested whether its expression is correlated with the expression of PARP1 in 498 primary neuroblastomas (GSE62564) and found a positive correlation between the two genes ( Figure 4B). Accordingly, E2F1 overexpression was also associated with advanced stage (Figure 4C), poor overall survival ( Figure 4D) and event-free survival ( Figure 4E) in the same dataset.
We also found a positive correlation between PARP1 and E2F1 protein levels in 6 selected NB cell lines. We observed that NB cells rs907187-CC have higher PARP1 (P = .0008) and E2F1 (P = .005) protein levels compared to rs907187-CG/GG cells ( Figure S6A-B).

| D ISCUSS I ON
The way in which a cancer patient responds to drug treatment is dependent on many variables and among these, the individual genotype is an essential factor for influencing drug behaviour.
Because of the importance of DNA repair process in determining drug sensitivity and resistance, in various cancer types, the role of polymorphisms in DNA repair genes have been explored to explain inter-individual differences in the treatment response or survival. 37  Here, we first show that high PARP1 mRNA levels associated with a worse outcome in three different datasets of patients with NB, suggesting that PARP1 expression may antagonize the effect of chemotherapeutic agents used in current NB therapy. We then selected and analysed two variants in PARP1 (rs907187 and rs2048426). Our data comprehensively reveal that the G allele of rs907187, located in promoter region, was associated with reduced activity of the dual-luciferase in a promoter reporter, with de-

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
MA performed the research and wrote the paper, VAL and FC analysed the data, SC, AC, AM and DC performed the research, RH, LA and AQ contributed essential reagents and tools, AG and AI designed the research study, MVC performed the research and analysed the data, and MC designed the research study and wrote the paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data generated and analysed during the current study are available from the corresponding author on reasonable request. Public data and data repositories are referenced within the manuscript.