Neuroblastoma (NB) is a pediatric tumor derived from the sympathoadrenal lineage of neural crest progenitor cells and represents the most common malignancy in early childhood . DNA and RNA aberrant profiles have been shown to identify mechanisms behind the clinical outcome of NB as the expression of several genes involved in proliferation, differentiation and metastasis that negatively impact on therapy success. Despite recent improvements in survival in randomized trials, nearly 50% of children with high-risk disease is refractory to therapy or suffer a relapse [2, 3, 4]. High-risk tumors are characterized by un-differentiated phenotype, age at diagnosis ≥18 months and harbor a very low rate of recurrent somatic mutations in both nuclear and mitochondrial DNA [5, 6, 7, 8, 9].
Hypoxia is an important factor in the pathology of many human diseases, including cancer, diabetes, aging, and stroke/ischemia. Low oxygen levels represent an important microenvironment condition to affect the activation status of signaling pathways as drug resistance mechanism. Indeed the increased expression of Hypoxia-Inducible-Factor HIF-1α mRNA (HIF1A) in tumors is relevant to establish resistance to therapeutic approaches as radiotherapy [10, 11]. We have recently reported that high HIF1A expression may stratify high-risk NB patients with poorer prognosis and low HIF1A expression enhances neuronal differentiation signaling pathways activation and response to differentiating agents . The identification of factors able to influence the expression levels of HIF1A could allow greater therapeutic success. Recent reports suggest that HIF1-α protein might be degraded in VHL-independent manner following intracellular accumulation of methylglioxal (MGO), a highly reactive α-oxoaldehyde formed as a by-product of glycolisis [13, 14]. Polymorphisms in glucoxylase I enzyme (GLOI) results in down-regulation of GLOI enzyme that play important role in MGO detoxyfication and favor damage from oxydative stress and the degradation pathway of HIF1A . Indeed, conditions with increased availability of glucose, such as diabetes or down-regulation of GLOI highlight the importance of mechanisms to disrupt cell response to hypoxia.
Tumor cells respond to repeated oxygen levels fluctuations in tumor microenvironment through epigenetic control. Epigenetic regulatory mechanisms are coordinated at several levels: i) DNA, by (hydroxy) methylation of CpG islands (CGI), ii) RNA, through involvement of regulatory noncoding RNA, and iii) proteins, by activation of epigenetic regulators and posttranslational modificators of histones. Their concerted action in hypoxia drives tumor plasticity through the acquisition of local or global chromatin modifications, which allow the accessibility of hypoxia-responsive elements (HRE) loci or of new active DNA regions at hypoxia inducible factors .
Epigenetic regulation of gene expression by DNA methylation plays a central role in determining tissue specific gene expression and chromosome instability. In cancer, the DNA methylation landscape is very complex: promoter CGIs hypermethylation is associated to inactivation of tumor suppressors as well as the presence of DNA hypomethylation blocks and contiguously hypermethylated CGIs at telomeric regions [17, 18]. Several studies show HIF1A expression can control DNA hypomethylation status of HRE. Interestingly, more than half of histone demethylase belonged to Jumonji C family genes were up-regulated by hypoxia and four of them (JMJD1A, JMJD2B, JMJD2C, PLU-1) were reported to be direct HIF1A targets and may result in increased HIF-1α binding to the HRE [19, 20].
Tumor hypoxia acts as a novel regulator of DNA methylation independently of HIF1A activity. High levels of hypoxia metabolites as succinate and fumarate altered the global DNA methylation patterns via significant DNA hypermethylation . Activity of ten-eleven translocation (TET) enzymes that catalyze DNA demethylation through 5-methylcytosine oxidation depends directly on oxygen shortage. Indeed, TETs activity is reduced by tumor hypoxia in human and mouse cells . Although HIF1A plays a role in defining DNA methylation status of its targets, its role in the global hypermethylation induced by hypoxia remains to be explored .
To shed light on the molecular mechanisms by which hypoxia reshapes gene expressions of tumors, we have performed an integrated analysis of gene expression and DNA methylation in NB cells upon HIF1A inhibition in normoxia and hypoxia conditions. We found that HIF1A transcription response in hypoxia is driven by epigenetic control of low oxygen levels and can upgrade high-risk tumor features. Interestingly, HIF1A targets expressed in both normoxic and hypoxic areas may provide novel targets to eradicate solid tumors.
The human SKNBE2 (ATCC #CRL-2271) cell line was grown in Dulbecco’s modified Eagle’s medium supplemented with 10% heat inactivated fetal bovine serum (Sigma), 1 mM L-glutamine, penicillin (100 U/ml) and streptomycin (100 μ g/ml) (Invitrogen), at 37 °C, under 5% CO2 in a humidified atmosphere. The cells exposed to hypoxia were grown at 0.5% oxygen for 2 h. The cells used for all the experiments were re-authenticated and tested as mycoplasma-free. Early-passage cells were used and cumulative culture length was less than 3 months after resuscitation.
Lentiviral production to knock-down HIF1A expression
To knock-down HIF1A expression, the pGIPZ lentiviral shRNAmir that targets human HIF1A were purchased from Open Biosystems (Thermo Fisher Scientific, Inc.). We used two different shRNAs against HIF1A: V2LHS_132152 (RHS4430–98513964) (shHIF1A#A) and V2LHS_236x718 (RHS443098513880) (shHIF1A#B). A non-silencing pGIPZ lentiviral shRNAmir was used as the control (RHS4346). The production of lentivirus particles and cells infection was performed as previously described . To obtain 100% GFP-positive cells, puromycin was added into the medium for an additional 10 days.
Fractionation of nuclear proteins and western blotting
Cell pellets were resuspended in a hypo-tonic buffer (10 mM HEPES-K +, pH 7.5, 10 mM KCl, 1.5 mM MgCl 2, 0.5 M dithiothreitol) in the presence of a protease inhibitors cocktail (Roche). The cells were lysed by addition of ice-cold 0.5% NP-40 for 10 min. The nuclei were pelleted at 1000 x g for 2 min at 4 °C and nuclear protein extraction and concentrations was determined as previously described . Protein membranes were probed with anti-HIF-1α (610,959; BD Biosciences) and horseradish-peroxidase-conjugated anti-mouse secondary antibody (1:4000 dilution; ImmunoReagent). Positive bands were visualized using the ECL kit SuperSignal West Pico Chemiluminescent Substrate (Pierce). A rabbit anti-H3 antibody (ab1791 Abcam) was used as the control for equal loading.
RNA isolation, cDNA library construction and sequencing
Total RNA was isolated from NB cell line using TRIzol LS Reagent (Invitrogen) according to manufacturer’s instructions; samples quality and library construction is described in Additional file 1. cDNA Sequencing was accomplished using an Illumina HiSeq™ 2000 platform according to the manufacturer’s protocols (Analysis performed at BIOGEM facility). Illumina paired end sequencing protocol yielded about 20 millions of 2x101nt reads with high quality bases (mean quality of 34) and mean % GC of 46.
Analysis of differentially expressed genes and gene set enrichment
Sequencing data were analyzed with the set of open source programs of the Tuxedo suite: TopHat v2.0.14 (for sequence alignment) and Cufflinks v2.1.0 (for differential expression analysis), following the pipeline published in Nature Protocols by Trapnell et al., 2012 (see Additional file 1 for further details) . The set of output files obtained by Cufflinks was inspected and explored using the R-Bioconductor package CummeRbund v2.16.0, which provides functions to read, subset, filter and plot results. We selected genes differentially expressed in each of the pairwise comparisons if the Benjamini Hochberg adjusted P-value (FDR) was under 0.05 and if the Log2 transformed fold change was greater than + 0.5 (up-regulated) or lower than − 0.5 (down-regulated). The lists of these genes were used to query Pathway and Gene Ontology databases. The functional enrichment analysis tool Webgestalt (WEB-based GEne SeT AnaLysis Toolkit) was used to detect significant enrichments for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The enrichment analysis was performed using the following criteria: an hypergeometric test for statistical analysis, FDR ≤ 0.05 and 10 as minimum number of genes for a category. Data generated during this study are included in this manuscript and in Additional files 1, 2 and 3.
DNA extraction, bisulfite modification and DNA methylation array hibridization
DNA extraction was carried out with the Wizard Genomic DNA Purification Kit (Promega, WI, USA), including a RNA removal step, according to the protocol provided by the supplier. The DNA was quantified with the Nanodrop and 1 μg was used for bisulfite modification using EZ-96 DNA Methylation™ Kit (Zymo Research CA, USA) with the modification step according to the recommendations for array processing of the samples. Control PCRs were carried out before array analysis to confirm successful modification of the DNA. The bisulfite-modified DNA (500 ng) was laid on the Infinium HumanMethylation450 BeadChips (Illumina), which determine the methylation levels of 485,000 CpG sites. The fluorescence signals were measured from the BeadArrays using an Illumina BeadStation GX scanner. The raw fluorescence images (IDAT files) were then analyzed using R and R-Bioconductor packages. The ChAMP package was used for data preprocessing, normalization and comparison between groups [25, 26]. Singular value decomposition analysis was performed to identify confounding factors and evaluate possible batch effects, while the SWAN (Subset Within Array Normalization) method was used for probe intensity normalization. After these steps the fraction of failed probes was about 0.003%.
Analysis of differentially methylated CpGs, CpGs enrichment and correlation to expression
ChAMP assigns a score called “β value” to each CpG site, which corresponds to the ratio between the fluorescence signal of the methylated allele (C) and the unmethylated (T) alleles. The β value, ranging from 0 to 1, represents the methylation status of each probe from totally unmethylated (β = 0) to totally methylated (β = 1). The software was used to calculate probes that were differentially methylated between groups. CpG sites were considered as differentially methylated, in a contrast, if Δβ (delta beta) was below − 0.2 (hypo-methylation) or above 0.2 (hyper-methylation) and the FDR was lower than 0.05. Data generated during this study are included in this manuscript and in Additional files 1, 2 and 3.
The expression levels of 13 genes were analyzed using real-time, quantitative PCR in SKNBE2 shHIF1A#A and shCTR cells. Total RNA extraction using TRIzol LS Reagent (Invitrogen) and cDNA retrotranscription using the High Capacity cDNA Reverse Transcription Script (Applied Biosistem) was performed according to the manufacturer protocol. The cDNA samples were diluted to 20 ng/μ l. Gene-specific primers were designed by using PRIMEREXPRESS software (Applied Biosystems) and primers sequences for each gene are listed in Additional file 1. Real-time PCR was performed using SYBR Green PCR Master Mix (AppliedBiosystems). All real-time PCR reactions were performed using the 7900HT Fast Real-Time PCR System (Applied Biosystems). The experiments were carried out in triplicate for each data point. The housekeeping gene β -actin was used as the internal control. Relative gene expression was calculated using the 2 −ΔΔCT method as described in our previous work , where the ∆CT was calculated using the differences in the mean CT between the selected genes and the internal control (β -actin). The mean fold change of 2 − (average ∆∆CT) was determined using the mean difference in the ∆CT between the gene of interest and the internal control.
HIF1A driven response in normoxia and in hypoxia conditions
SKNBE2 NB cells have biochemical features of neurons and display NMYC amplification, a marker of NB malignant progression. SKNBE2 cells were depleted for HIF1A expression (shHIF1A) by the use of two short hairpin against HIF1A (SKNBE2 shHIF1A#A and SKNBE2 shHIF1A#B) and were grown in normoxia and hypoxia conditions (NX and HYP); unsilenced cells were used as control (SKNBE2 shCTR) (Fig. 1a). To evaluate the hypoxic status of the cells after their exposure to low oxygen conditions we tested the expression of known hypoxia targets (Additional file 2: Figure S1). To provide genes and pathways differentially regulated by HIF1A triplicates of silenced (SKNBE2 shHIF1A#B NX and SKNBE2 shHIF1A#B HYP) and unsilenced (SKNBE2 shCTR NX and SKNBE2 shCTR HYP) cells were subjected to RNA-seq experiment.
To get insights into the relationships between the experimental conditions, we performed hierarchical clustering of gene-based FPKM counts. Results clearly showed two main branches of the dendrogram separating silenced and unsilenced conditions (Additional file 2: Figure S2A). Although small changes of global expression levels were observed in the four conditions (Additional file 2: Figure S2B), principal component analysis (PCA) and multi-dimensional scaling (MDS) highlighted the important role of oxygen in separating shCTR HYP and shCTR NX whereas shHIF1A NX and shHIF1A HYP showed highly similar profiles in PCA analysis (Additional file 2: Figure S2C and D).
By comparing gene expression in shHIF1A NX vs shCTR NX and in shHIF1A HYP vs shCTR HYP, we obtained two gene sets. Raw calls of differentially expressed genes were subsequently filtered by fold change (Log2 ≥ + 0.5 or ≤ − 0.5) and statistical significance (FDR ≤ 0.05) (Additional file 3: Tables S1 and S2). HIF1A silencing in normoxia affects the expression of much more genes (“shHIF1A NX vs shCTR NX” includes 2656 genes) than in hypoxia (“shHIF1A HYP vs shCTR HYP” includes 1886 genes) (Additional file 2: Figure S3A). KEGG pathway analysis showed that the most significantly enriched terms were “metabolic pathway” in shHIF1A NX vs shCTR NX and “axon guidance” in shHIF1A HYP vs shCTR HYP (FDR ≤ 0.05) (Additional file 2: Figure S3B).
By intersecting the above-cited “two gene sets”, we obtained three gene lists: 1) genes (n = 630) regulated “exclusively in hypoxia”; 2) genes (n = 1400) regulated “exclusively in normoxia” and 3) a list of genes (n = 1256) which are commonly regulated by HIF1A that we named HIF1A target genes (Fig. 1b). The expression trend of 1237 out of 1256 HIF1A target genes (98.88%) was concordant upon HIF1A depletion in NX and HYP (Fig. 1c). Conversely, 19 genes out of 1256 (less than 1.2%) had an opposite regulation in shHIF1A NX vs shCTR NX and shHIF1A HYP vs shCTR HYP, suggesting they are downstream targets of HIF1A related pathways. KEGG pathway analysis of the two “exclusive” gene sets revealed an enrichment of metabolic pathways in normoxia and an enrichment of axon guidance and pathways in cancer in hypoxia (FDR ≤ 0.05). KEGG pathway analysis of HIF1A target genes revealed an enrichment of MAPK signaling pathway, pathways in cancer and axon guidance (FDR ≤ 0.05) (Fig. 1d).
The reliability of RNA-seq data was assessed by RT-PCR in SKNBE2 shCTR and shHIF1A#A (Fig. 1e). We validated genes that have RNAseq log2 fold change ranging from − 2 to 2, in shCTR HYP vs shCTR NX (Additional file 3: Table S3) and shHIF1A HYP vs shCTR HYP gene list. We found that expression levels measured by RNA-Seq were consistent with those obtained by RT-PCR.
Additionally, we confirmed these results by RT-PCR in SHSY5Y NB cells that have biochemical features of neurons but do not display high-risk marker as NMYC amplification. As described in Supplementary data, SHSY5Y cells were depleted for HIF1A expression (shHIF1A) and unsilenced cells were used as control (shCTR). The gene expression levels measured by RT-PCR in SHSY5Y cells were consistent with those obtained by RNA-Seq in SKNBE2 cells (Additional file 2: Figure S4).
Transcription activity under hypoxia exposure is both HIF1A dependent and HIF1A independent
The above results clearly highlight that HIF1A has a diverse role in normoxia than in hypoxia. Our hypothesis is that HIF1A driven response depends on the epigenetic reprogramming caused by low oxygen levels that may shape chromatin state and give HIF1A accessibility to HRE DNA regions previously closed. Furthermore, chromatin may remodel in regions not comprising HIF1A targets and allow HIF1A and other transcription factors access to new active DNA regions. To deeply investigate which genes are HIF1A dependently and HIF1A independently expressed, gene set shHIF1A HYP vs shCTR HYP (n = 1886 gene, Additional file 3: Table S2) and gene set of shCTR HYP vs shCTR NX (n = 3263 genes, Additional file 3: Table S3) were crossed. We found that 674 genes were regulated in both gene sets (Fig. 2a). Interestingly, 420 out of 674 genes show the same regulation in both gene sets (Log2); the expression of these genes named “Hypoxia targets” is affected by low oxygen concentrations and not affected by HIF1A (Fig. 2b). KEGG pathway analysis shows that “Hypoxia targets” are enriched in pathway that regulate cytoskeleton, ligand-receptor interaction and axon guidance (FDR ≤ 0.05). By contrast, 254 out of 674 genes have an opposite regulation in the two lists (Log2). These genes might be direct targets of HIF1A (here named: “HIF1A direct-targets”) in hypoxia because when HIF1A is depleted their regulation is inverted (Fig. 2c). KEGG pathway analysis shows that “HIF1A direct-targets” are enriched in metabolic and cancer pathways similar to pathways affected by HIF1A silencing “exclusively in hypoxia” (FDR ≤ 0.05). These findings suggest that NB cells adapt to hypoxia by HIF1A-dependent and HIF1A-independent driven response.
DNA sites have variable methylation status under hypoxia exposure
Genome-wide methylation analysis using Infinium HumanMethylation450 BeadChips was performed in triplicate as described for RNAseq. To get an overview of the methylation patterns in the normalized data, hierarchical clustering of the most variable probes was performed. The analysis separated samples into four clusters, one for each experimental condition, within which replicates are grouped (Fig. 3a). Probes hypo or hyper-methylated with a Δβ-value greater than 0.2 (20%) in shHIF1A HYP vs shCTR HYP (named HIF1A probes) and in shCTR HYP vs shCTR NX (named Hypoxia probes) were selected. The sets of HIF1A probes and Hypoxia probes include 1078 (Additional file 3: Table S4) and 260 (Additional file 3: Table S5) differentially methylated CpG sites respectively. A global hypermethylation status of Hypoxia probes (Δβ ≥ 0.2) and hypomethylation status of HIF1A probes (Δβ ≤ − 0.2) was observed (Fig. 3b). Both probe sets cluster close to each other (Fig. 3a) whereas probes differentially methylated in shHIF1A NX vs shCTR NX were not observed.
We measured the genomic distribution of Hypoxia probes and HIF1A probes in relation to CGI centric annotation. Hypoxia probes were overrepresented in regions of low CpGs (open sea, > 4Kb from the CGI) and underrepresented in CGIs (island). Once HIF1A is depleted (HIF1A probes) open sea and island regions became under- and over-represented, respectively, (Fig. 3c). CGI shelves (>2Kb from the CGI) and shores (<2Kb from the CGI) showed a non-random distribution between the two probe sets. Hypoxia probes and HIF1A probes were mapped in relation to gene positions. In both contrasts, a strong overrepresentation of probes in intergenic regions (IGR) and gene body as well as an overall underrepresentation of probes located in the first exon, 3′ and 5’UTR and in the upstream region of transcription site (TSS) was observed (Fig. 3d). Hypoxia probes and HIF1A probes were crossed and 150 probes were found methylated in both probes lists. The methylation status (Δβ) of 150 common probes in hypoxia (shCTR HYP vs shCTR NX, Hypoxia probes) is reverted once HIF1A is depleted (shHIF1A HYP vs shCTR HYP, HIF1A probes) (Fig. 3e) and a strong overrepresentation of common probes in IGR is observed (Fig. 3f). Overall, these results suggest DNA methylation status is strictly correlated with oxygen storage and HIF1A control of DNA methylation of IGR, gene body and TSS probes could occur only upon hypoxia induced epigenetic reprogramming.
Correlation of DNA methylation and gene expression under hypoxia exposure
The correlation between the differential expression and differential methylation was explored for each gene-probe pair (p-value ≤0.05). We searched for changes in opposite directions (eg. Up-regulation of the gene expression and hypo-methylation of the related CpG probe). In the shHIF1A HYP vs shCTR HYP comparison, we selected 31 gene-probe pairs (p ≤ 0.05, Fig. 4a, Table 1), whereas in the shCTR HYP vs shCTR NX comparison 18 gene-probe pairs were kept (p ≤ 0.05, Fig. 4b, Table 1). These probes are located on regulatory regions as UTR, TSS200, TSS1500 and in gene bodies (Additional file 2: Figure S5).
Correlation between DNA methylation (Probe ID) and gene expression (Gene) in hypoxia
|ProbeID||Gene||Δβ HIF1A||Log2 FC HIF1A||Gene feature||CpG feature|
|ProbeID||Gene||Δβ Hypoxia||Log2 FC Hypoxia||Gene feature||CpG feature|
Δβ Hypoxia: Delta beta in shCTR HYP vs shCTR NX.; Δβ HIF1A: Delta beta in shHIF1A HYP vs shCTR HYP; Log2 FC: log2 of expression fold change
We explored gene-probe pairs correlation in 105 NB tumors for which matched methylation and gene expression data were available (GEO accessions: GSE73515 and GSE73517, respectively) and restricted our analysis to Low risk (n = 40) and High risk (n = 56) tumors as defined by Henrich et al. . We found that the correlations for KIF26B, EFCAB2, BCL2L11, VAV2 and SORBS2 gene expression with methylation status were validated (Additional file 2: Figure S6A-E; P < 0.05). Additionally, in an independent set of NB tumors (GSE16476), we found that the expression of these genes was also associated to NB patient’s survival (Additional file 2: Figure S6F).
The two gene signatures generated from gene-probe pairs were named “HIF1A signature in hypoxia” and “Hypoxia signature”. To assess the prognostic potential of these signatures, we used Low risk (n = 40) and High risk (n = 56) tumors from the GSE73517 gene expression data set . Hierarchical clustering based on Euclidean distances of expression levels, showed that the genes in “HIF1A signature in hypoxia” (Fig. 4c) clustered the 35.7% (20/56) of High risk and the 77.5% (31/40) of Low risk patients, in two separate groups (P < 1.0 × 10− 4) according to their Risk category. In contrast, “Hypoxia signature” (Fig. 4d), clustered the 50% (28/56) of High risk and the 57.5% (23/40) of Low risk patients according to their Risk category (P < 1.0 × 10− 4). Indeed, we verified that our gene signatures correctly classified a discrete portion of both High and Low risk patients.
Identification of enhancers methylated under low oxygen conditions
The strong enrichment of hypoxia differentially methylated sites in IGR suggests that the changes of methylation pattern mainly occur at putative regulatory regions distant from target genes (as shown in Fig. 3d-f). To further select which hypoxia differentially methylated sites (Hypoxia probes n = 260) are located in NB putative regulatory regions we re-analysed DNase hypersensitivity assay and Chip-Seq histone acetylation (H3K27ac) data deposited in Gene Expression Omnibus database (GSE65664, Additional file 3: Table S6) for additional SKNBE2, CHP134 and SHSY5Y cell lines. We obtained that 113 out of 260 probes (43.5%) were located in regulatory active regions and 14 probes were annotated in all cell lines with at least 4 epigenetic markers (Additional file 3: Table S6). We searched the genes distant 1 Mb up- or down-stream from the 14 aforementioned probes in the RNAseq data (shCTR HYP vs shCTR NX, Log2 ≥ 0.2, ≤ − 0.2) and we listed the candidate targets of these regulative regions in Table 2 (and Additional file 3: Table S7). We observed that the methylation status of 9 out of 14 probes is affected by oxygen levels (Δβ Hypoxia; Hypoxia probes) but not by HIF1A expression (Δβ HIF1A; HIF1A probes). Putative targets show gene expression levels affected by oxygen levels (Log FC Hypoxia, in shCTR HYP vs shCTR NX, Log2 ≥ 0.2, ≤ − 0.2) and not by HIF1A expression (Log FC HIF1A, in shHIF1A HYP vs shCTR HYP, Log2 ≥ 0.2, ≤ − 0.2). Conversely, the methylation status of the remaining 5 probes is inversely affected in both conditions (Δβ Hypoxia and Δβ HIF1A) as the expression of putative targets (Log FC Hypoxia and Log FC HIF1A). The best NB candidate targets are represented by genes (highlighted in Table 2) whose expression correlates with NB patient’s event-free survival (Additional file 2: Figure S7).
Increased expression of HIF1A in tumors is relevant to establish resistance to therapy [10, 11]. Interestingly, we have previously reported that high HIF1A expression may stratify high-risk NB patients with poorer prognosis .
Currently, targeting of hypoxia signaling has limitations in clinics with regard to changeable oxygen concentrations in solid tumor areas and HIF1A direct compounds do not show clinical efficiency. Indeed, the identification of HIF1A target genes and deep insights into the mechanisms of HIF1A driven gene expression may provide novel risk factors to meliorate survival/therapeutic successes in patients with high-risk tumors that lack of precisely genomic causes.
In the present study, we have investigated HIF-1 driven transcription activity in both hypoxic and normoxic conditions in NB cells depleted of HIF1A expression. The analysis of pathways regulated by HIF1A exclusively in normoxic NB cells shows a role of HIF1A in metabolic process necessary for tumor cells viability. Particularly, the global down-regulation of gene expression in absence of HIF1A suggests that NB cells slow down their metabolic activity, thus becoming less proliferating. HIF1A involvement in basic cellular activity, like glycolytic pathways, has been described .
Contrary, in hypoxic cells the absence of HIF1A affects the activation of neuronal differentiation pathways in line with literature data showing that low oxygen in environments causes de-differentiation of NB cells towards an immature and neural-crest-like phenotype . We have previously highlighted HIF1A involvement in NB neuronal differentiation pathways activation and response to differentiating agents .
Interesting to note, mostly of genes regulated by HIF1A in both normoxic and hypoxic areas belong to MAPK pathways. This pathway is frequently altered in high-risk NB at relapse and at diagnosis and multiple drugs aimed to target MAPK signaling are used in current clinical trials for the treatment of metastatic tumors [5, 8, 31]. Indeed, HIF1A target genes in both normoxic and hypoxic areas may provide potential targets for a precision therapy. HIF1A is not the unique player to define the whole picture of hypoxia-regulated gene expression. In effect, we report that NB cells adapt to hypoxia by HIF1A-dependend and HIF1A-indipendent driven response. These findings help us to understand how oxygen is sensed at NB cellular levels.
We assume that HIF1A driven transcriptional response in hypoxia is a consequence of the epigenetic control of low oxygen levels at DNA methylation status. We have observed that hypoxia exposure induces a global DNA hypermethylation in NB cells and HIF1A itself might control DNA methylation status. A global DNA hypermethylation has been previously linked to poor NB prognosis as site-specific DNA hypermethylation of tumor suppressor genes to optimize the environment for cancer initiation and progression [32, 33]. The hypoxia epigenetic controls at the levels of RNA and proteins still remain to be explored.
Despite the stereotype, DNA methylation does not appear to play a major role in gene regulation from 5’CGI promoters of most genes in hypoxia. Indeed, few genes show a correlation between expression and methylation status of close regulatory regions and some correlations were validated in NB samples. Hypoxic gene signatures generated from this correlation analysis are able to stratify NB patients in two risk categories. Although numerous prognostic gene signatures have been developed to classify NB patients, none has been introduced into clinical risk stratification systems [2, 34, 35]. To overcome these limitations, the establishment of gene signatures that take into account the effects of oxygen levels in tumor bulk more than clinical and genetic markers may be an innovative strategy for NB stratification at diagnosis. Of course, these findings need independent validations.
Conversely, low oxygen levels and HIF1A affect the methylation status of probes located in intragenic and intergenic regions [36, 37, 38]. Most probes are located in NB active regulatory regions and the different methylation status correlates to different expression of distant candidate targets associated with NB survival. These genes have been previously associated to therapy resistance and cancer progression and may represent potential markers for NB.
CDC20 is a component of the mammalian cell-cycle mechanism and activates the anaphase-promoting complex (APC); its inhibition may enhance radio sensitivity in nasopharyngeal carcinoma cells .
SNRPE (small nuclear ribonucleoprotein polypeptide E) has oncogenic effects in prostate cancer . TDP1 (Tyrosyl-DNA Phosphodiesterase 1) is DNA repair enzyme potential therapeutic target for the treatment of colorectal cancer .
FOXM1 (Forkhead Box M1) transcription factor regulates the expression of cell cycle genes and plays an important role in NB tumorigenicity through maintenance of cells undifferentiated state . Interestingly, FOXM1 overexpression in hypoxia has been already documented in cancer .
DMAP1 (DNA Methyltrasferase 1 Associated Protein 1) contributes to epatocarcinoma malignancy .
YBEY (C21orf57) is a highly conserved metalloprotein not-well characterized in cancer.
High-throughput sequencing-based studies have shown low mutations frequency in coding-portion of NB genome and high recurrence of structural rearrangement. Previous genome-wide association studies revealed that many loci associated with NB susceptibility lie in non-coding regions of the genome [35, 44, 45, 46]. Based on these evidences, it is reasonable to expect that recurrent non-coding somatic mutations could have a regulatory effect in NB tumorigenesis. In light of all this, our results further underline the role of non-coding regulatory elements in driving NB tumorigenesis through epigenetic regulation in hypoxia. How epigenetic landscape in hypoxia contributes to transformations and how these alterations complement other acquired somatic mutations need to be elucidated.
One limitation of this study is the use of established cell lines that reflects limited aspects of in vivo tumor microenvironments. It lacks geometrical complexity, cellular components including immune cells and organ-specific stromal cells, and extracellular matrix components. Here, our aim was to establish a HIF1A-based method useful in the investigation of undiscovered mechanisms of neuroblastoma tumorigenesis under hypoxic microenvironments. However, our results still need to be confirmed by functional validation and mechanistic studies which could further improve in vitro cell line models predictive validity.
Recent evidences of HIF1A gene expression increment in solid tumors and the stabilization of low HIF1A protein levels in normoxic cells highlight HIF1A transcription activity in both normoxic and hypoxic conditions. In the present study, we have investigated HIF1A targets activated in both oxygen level conditions by an analysis of gene expression and DNA methylation of NB cells silenced or unsilenced for HIF1A expression. We have verified that HIF1A transcription activity depends on oxygen levels and HIF1A targets regulated in both conditions might provide potential therapeutic targets to eradicate solid tumors. Hypoxia signatures might provide novel risk factors for NB stratification at diagnosis. Hypoxia regulates gene expression through an epigenetic control on regulatory elements distant from target genes (Fig. 5). Overall, the presented results may help to understand the molecular mechanisms by which hypoxia reshapes tumors and provide new direction for hypoxia solid tumor treatment.
HIF-1: Heterodimeric transcription factor hypoxia-inducible factor 1
HIF1A: Hypoxia-inducible factor α
We thank Fondazione Umberto Veronesi and AIRC.
This study was supported by grants from Associazione Italiana per la Ricerca sul Cancro (AIRC) (20757 and 19255); Ministero della Salute (GR-2011-02348722); “Fondazione Italiana per la Lotta al Neuroblastoma”; OPEN Associazione Oncologica Pediatrica e Neuroblastoma. : Regione Campania “SATIN” grant 2018-2020. F.C. was supported by Fondazione Umberto Veronesi post-Doc Fellowship. In all cases, the funders had no role in study design, data collection, analysis and interpretation of data, or in the writing of the manuscript and decision to submit it for publication.
Availability of data and materials
Data generated during this study are included in this manuscript and in supplementary information files.
FC conceived the idea, analysed and interpreted data, and drafted the manuscript. FC is responsible for the integrity of the work as a whole. MA and LP made substantial contributions to the RNA and DNA sample preparation, RTP-PCR and western blotting. VL performed bioinformatic analysis of RNAseq and DNA methylation data. AM, FV and MVC performed validation analysis. MC and AI edited the manuscript. All authors have read and approved the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Publisher’s Note (BMC)
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access – HIF-1 transcription activity….
Additional file 1: Supplementary Material and Methods cDNA library construction; Analysis of differentially expressed genes; Correlation of DNA methylation and gene expression; Primers sequences for RT-PCR. (DOC 55 kb)Additional file 2:Figure S1. Main features of RNA-Seq data; Figure S2.Differentially expressed genes in shHIF1A NX vs shCTR NX and in shHIF1A HYP vs shCTR HYP gene sets; Figure S3. Validation of RNAseq data by RT-PCR in SHSY5Y cells; Figure S4. Correlation of DNA methylation and gene expression; Figure S5.Genes differentially expressed and methylated under hypoxia and associated with NB survival; Figure S6. Survival analysis of NB patients based on the gene expression of candidate targets of putative enhancers differentially regulated in hypoxia. Figure S7. Survival analysis of NB patients based on the gene expression of candidate targets of putative enhancers differentially regulated in hypoxia. (DOCX 7156 kb)Additional file 3:Table S1. shHIF1A NX vs shCTR NX; Table S2. shHIF1A HYP vs shCTR HYP; Table S3. shCTR HYP vs shCTR NX; Table S4. Probes shCTR HYP vs shHIF1A HYP; Table S5. Probes shCTR HYP vs shCTR NX; Table S6. NB annotated enhancers; Table S7. Genes surrounding 14 probes. (XLS 1732 kb)
- Maris JM. Recent advances in neuroblastoma. N Engl J Med. 2010;362:2202–11.View ArticleGoogle Scholar
- Formicola D, Petrosino G, Lasorsa VA, Pignataro P, Cimmino F, Vetrella S, Longo L, Tonini GP, Oberthuer A, Iolascon A, Fischer M, Capasso M. An 18 gene expression-based score classifier predicts the clinical outcome in stage 4 neuroblastoma. J Transl Med. 2016;14:142.View ArticleGoogle Scholar
- Jung M, Russell AJ, Liu B, George J, Liu PY, Liu T, DeFazio A, Bowtell DD, Oberthuer A, London WB, Fletcher JI, Haber M, Norris MD, Henderson MJ. A Myc activity signature predicts poor clinical outcomes in Myc-associated cancers. Cancer Res. 2017;77:971–81.View ArticleGoogle Scholar
- Schleiermacher G, Mosseri V, London WB, Maris JM, Brodeur GM, Attiyeh E, Haber M, Khan J, Nakagawara A, Speleman F, Noguera R, Tonini GP, Fischer M, Ambros I, Monclair T, Matthay KK, Ambros P, Cohn SL, Pearson AD. Segmental chromosomal alterations have prognostic impact in neuroblastoma: a report from the INRG project. Br J Cancer. 2012;107:1418–22.View ArticleGoogle Scholar
- Pugh TJ, Morozova O, Attiyeh EF, Asgharzadeh S, Wei JS, Auclair D, Carter SL, Cibulskis K, Hanna M, Kiezun A, Kim J, Lawrence MS, Lichenstein L, McKenna A, Pedamallu CS, Ramos AH, Shefler E, Sivachenko A, Sougnez C, Stewart C, Ally A, Birol I, Chiu R, Corbett RD, Hirst M, Jackman SD, Kamoh B, Khodabakshi AH, Krzywinski M, Lo A, Moore RA, Mungall KL, Qian J, Tam A, Thiessen N, Zhao Y, Cole KA, Diamond M, Diskin SJ, Mosse YP, Wood AC, Ji L, Sposto R, Badgett T, London WB, Moyer Y, Gastier-Foster JM, Smith MA, Guidry Auvil JM, Gerhard DS, Hogarty MD, Jones SJ, Lander ES, Gabriel SB, Getz G, Seeger RC, Khan J, Marra MA, Meyerson M, Maris JM. The genetic landscape of high-risk neuroblastoma. Nat Genet. 2013;45:279–84.View ArticleGoogle Scholar
- Calabrese FM, Clima R, Pignataro P, Lasorsa VA, Hogarty MD, Castellano A, Conte M, Tonini GP, Iolascon A, Gasparre G, Capasso M. A comprehensive characterization of rare mitochondrial DNA variants in neuroblastoma. Oncotarget. 2016;7:49246–58.View ArticleGoogle Scholar
- Valentijn LJ, Koster J, Zwijnenburg DA, Hasselt NE, van Sluis P, Volckmann R, van Noesel MM, George RE, Tytgat GA, Molenaar JJ, Versteeg R. TERT rearrangements are frequent in neuroblastoma and identify aggressive tumors. Nat Genet. 2015;47:1411–4.View ArticleGoogle Scholar
- Eleveld TF, Oldridge DA, Bernard V, Koster J, Colmet Daage L, Diskin SJ, Schild L, Bentahar NB, Bellini A, Chicard M, Lapouble E, Combaret V, Legoix-Ne P, Michon J, Pugh TJ, Hart LS, Rader J, Attiyeh EF, Wei JS, Zhang S, Naranjo A, Gastier-Foster JM, Hogarty MD, Asgharzadeh S, Smith MA, Guidry Auvil JM, Watkins TB, Zwijnenburg DA, Ebus ME, van Sluis P, Hakkert A, van Wezel E, van der Schoot CE, Westerhout EM, Schulte JH, Tytgat GA, Dolman ME, Janoueix-Lerosey I, Gerhard DS, Caron HN, Delattre O, Khan J, Versteeg R, Schleiermacher G, Molenaar JJ, Maris JM. Relapsed neuroblastomas show frequent RAS-MAPK pathway mutations. Nat Genet. 2015;47:864–71.View ArticleGoogle Scholar
- Lasorsa VA, Formicola D, Pignataro P, Cimmino F, Calabrese FM, Mora J, Esposito MR, Pantile M, Zanon C, De Mariano M, Longo L, Hogarty MD, de Torres C, Tonini GP, Iolascon A, Capasso M. Exome and deep sequencing of clinically aggressive neuroblastoma reveal somatic mutations that affect key pathways involved in cancer progression. Oncotarget. 2016;7:21840–52.View ArticleGoogle Scholar
- Jia X, Hong Q, Lei L, Li D, Li J, Mo M, Wang Y, Shao Z, Shen Z, Cheng J, Liu G. Basal and therapy-driven hypoxia-inducible factor-1alpha confers resistance to endocrine therapy in estrogen receptor-positive breast cancer. Oncotarget. 2015;6:8648–62.PubMedPubMed CentralGoogle Scholar
- Jiang YZ, Liu YR, Xu XE, Jin X, Hu X, Yu KD, Shao ZM. Transcriptome analysis of triple-negative breast Cancer reveals an integrated mRNA-lncRNA signature with predictive and prognostic value. Cancer Res. 2016;76:2105–14.View ArticleGoogle Scholar
- Cimmino F, Pezone L, Avitabile M, Acierno G, Andolfo I, Capasso M, Iolascon A. Inhibition of hypoxia inducible factors combined with all-trans retinoic acid treatment enhances glial transdifferentiation of neuroblastoma cells. Sci Rep. 2015;5:11158.View ArticleGoogle Scholar
- Inagi R, Kumagai T, Fujita T, Nangaku M. The role of glyoxalase system in renal hypoxia. Adv Exp Med Biol. 2010;662:49–55.View ArticleGoogle Scholar
- Bento CF, Fernandes R, Ramalho J, Marques C, Shang F, Taylor A, Pereira P. The chaperone-dependent ubiquitin ligase CHIP targets HIF-1alpha for degradation in the presence of methylglyoxal. PLoS One. 2010;5:e15062.View ArticleGoogle Scholar
- Donato L, Scimone C, Nicocia G, Denaro L, Robledo R, Sidoti A, D’Angelo R. GLO1 gene polymorphisms and their association with retinitis pigmentosa: a case-control study in a Sicilian population. Mol Biol Rep. 2018;45:1349–55.View ArticleGoogle Scholar
- Watson JA, Watson CJ, McCann A, Baugh J. Epigenetics, the epicenter of the hypoxic response. Epigenetics. 2010;5:293–6.View ArticleGoogle Scholar
- Buckley PG, Das S, Bryan K, Watters KM, Alcock L, Koster J, Versteeg R, Stallings RL. Genome-wide DNA methylation analysis of neuroblastic tumors reveals clinically relevant epigenetic events and large-scale epigenomic alterations localized to telomeric regions. Int J Cancer. 2011;128:2296–305.View ArticleGoogle Scholar
- Doi A, Park IH, Wen B, Murakami P, Aryee MJ, Irizarry R, Herb B, Ladd-Acosta C, Rho J, Loewer S, Miller J, Schlaeger T, Daley GQ, Feinberg AP. Differential methylation of tissue- and cancer-specific CpG island shores distinguishes human induced pluripotent stem cells, embryonic stem cells and fibroblasts. Nat Genet. 2009;41:1350–3.View ArticleGoogle Scholar
- Koslowski M, Luxemburger U, Tureci O, Sahin U. Tumor-associated CpG demethylation augments hypoxia-induced effects by positive autoregulation of HIF-1alpha. Oncogene. 2011;30:876–82.View ArticleGoogle Scholar
- Xia X, Lemieux ME, Li W, Carroll JS, Brown M, Liu XS, Kung AL. Integrative analysis of HIF binding and transactivation reveals its role in maintaining histone methylation homeostasis. Proc Natl Acad Sci U S A. 2009;106:4260–5.View ArticleGoogle Scholar
- Wentzel JF, Lewies A, Bronkhorst AJ, van Dyk E, du Plessis LH, Pretorius PJ. Exposure to high levels of fumarate and succinate leads to apoptotic cytotoxicity and altered global DNA methylation profiles in vitro. Biochimie. 2017;135:28–34.View ArticleGoogle Scholar
- Thienpont B, Steinbacher J, Zhao H, D’Anna F, Kuchnio A, Ploumakis A, Ghesquiere B, Van Dyck L, Boeckx B, Schoonjans L, Hermans E, Amant F, Kristensen VN, Peng Koh K, Mazzone M, Coleman M, Carell T, Carmeliet P, Lambrechts D. Tumour hypoxia causes DNA hypermethylation by reducing TET activity. Nature. 2016;537:63–8.View ArticleGoogle Scholar
- Wu CY, Tsai YP, Wu MZ, Teng SC, Wu KJ. Epigenetic reprogramming and post-transcriptional regulation during the epithelial-mesenchymal transition. Trends Genet. 2012;28:454–63.View ArticleGoogle Scholar
- Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and cufflinks. Nat Protoc. 2012;7:562–78.View ArticleGoogle Scholar
- Butcher LM, Beck S. Probe lasso: a novel method to rope in differentially methylated regions with 450K DNA methylation data. Methods. 2015;72:21–8.View ArticleGoogle Scholar
- Morris TJ, Butcher LM, Feber A, Teschendorff AE, Chakravarthy AR, Wojdacz TK, Beck S. ChAMP: 450k Chip analysis methylation pipeline. Bioinformatics. 2014;30:428–30.View ArticleGoogle Scholar
- Capasso M, Avvisati RA, Piscopo C, Laforgia N, Raimondi F, de Angelis F, Iolascon A. Cytokine gene polymorphisms in Italian preterm infants: association between interleukin-10 -1082 G/a polymorphism and respiratory distress syndrome. Pediatr Res. 2007;61:313–7.View ArticleGoogle Scholar
- Henrich KO, Bender S, Saadati M, Dreidax D, Gartlgruber M, Shao C, Herrmann C, Wiesenfarth M, Parzonka M, Wehrmann L, Fischer M, Duffy DJ, Bell E, Torkov A, Schmezer P, Plass C, Hofer T, Benner A, Pfister SM, Westermann F. Integrative genome-scale analysis identifies epigenetic mechanisms of transcriptional deregulation in unfavorable neuroblastomas. Cancer Res. 2016;76:5523–37.View ArticleGoogle Scholar
- Hu CJ, Iyer S, Sataur A, Covello KL, Chodosh LA, Simon MC. Differential regulation of the transcriptional activities of hypoxia-inducible factor 1 alpha (HIF-1alpha) and HIF-2alpha in stem cells. Mol Cell Biol. 2006;26:3514–26.View ArticleGoogle Scholar
- Das B, Tsuchida R, Malkin D, Koren G, Baruchel S, Yeger H. Hypoxia enhances tumor stemness by increasing the invasive and tumorigenic side population fraction. Stem Cells. 2008;26:1818–30.View ArticleGoogle Scholar
- Bockorny B, Rusan M, Chen W, Liao RG, Li Y, Piccioni F, Wang J, Tan L, Thorner AR, Li T, Zhang Y, Miao C, Ovesen T, Shapiro GI, Kwiatkowski DJ, Gray NS, Meyerson M, Hammerman PS, Bass AJ. RAS-MAPK reactivation facilitates acquired resistance in FGFR1-amplified lung Cancer and underlies a rationale for upfront FGFR-MEK blockade. Mol Cancer Ther. 2018;17:1526–39.View ArticleGoogle Scholar
- Olsson M, Beck S, Kogner P, Martinsson T, Caren H. Genome-wide methylation profiling identifies novel methylated genes in neuroblastoma tumors. Epigenetics. 2016;11:74–84.View ArticleGoogle Scholar
- Sandoval J, Esteller M. Cancer epigenomics: beyond genomics. Curr Opin Genet Dev. 2012;22:50–5.View ArticleGoogle Scholar
- Asgharzadeh S, Pique-Regi R, Sposto R, Wang H, Yang Y, Shimada H, Matthay K, Buckley J, Ortega A, Seeger RC. Prognostic significance of gene expression profiles of metastatic neuroblastomas lacking MYCN gene amplification. J Natl Cancer Inst. 2006;98:1193–203.View ArticleGoogle Scholar
- Capasso M, Diskin SJ. Genetics and genomics of neuroblastoma. Cancer Treat Res. 2010;155:65–84.View ArticleGoogle Scholar
- Ehrlich KC, Paterson HL, Lacey M, Ehrlich M. DNA Hypomethylation in intragenic and intergenic enhancer chromatin of muscle-specific genes usually correlates with their expression. Yale J Biol Med. 2016;89:441–55.PubMedPubMed CentralGoogle Scholar
- Maunakea AK, Nagarajan RP, Bilenky M, Ballinger TJ, D’Souza C, Fouse SD, Johnson BE, Hong C, Nielsen C, Zhao Y, Turecki G, Delaney A, Varhol R, Thiessen N, Shchors K, Heine VM, Rowitch DH, Xing X, Fiore C, Schillebeeckx M, Jones SJ, Haussler D, Marra MA, Hirst M, Wang T, Costello JF. Conserved role of intragenic DNA methylation in regulating alternative promoters. Nature. 2010;466:253–7.View ArticleGoogle Scholar
- Neri F, Rapelli S, Krepelova A, Incarnato D, Parlato C, Basile G, Maldotti M, Anselmi F, Oliviero S. Intragenic DNA methylation prevents spurious transcription initiation. Nature. 2017;543:72–7.View ArticleGoogle Scholar
- Wang C, Su Z, Hou H, Li D, Pan Z, Tian W, Mo C. Inhibition of anaphase-promoting complex by silence APC/C (Cdh1) to enhance Radiosensitivity of nasopharyngeal carcinoma cells. J Cell Biochem. 2017;118:3150–7.View ArticleGoogle Scholar
- Anchi T, Tamura K, Furihata M, Satake H, Sakoda H, Kawada C, Kamei M, Shimamoto T, Fukuhara H, Fukata S, Ashida S, Karashima T, Yamasaki I, Yasuda M, Kamada M, Inoue K, Shuin T. SNRPE is involved in cell proliferation and progression of high-grade prostate cancer through the regulation of androgen receptor expression. Oncol Lett. 2012;3:264–8.View ArticleGoogle Scholar
- Meisenberg C, Gilbert DC, Chalmers A, Haley V, Gollins S, Ward SE, El-Khamisy SF. Clinical and cellular roles for TDP1 and TOP1 in modulating colorectal cancer response to irinotecan. Mol Cancer Ther. 2015;14:575–85.View ArticleGoogle Scholar
- Wang Z, Park HJ, Carr JR, Chen YJ, Zheng Y, Li J, Tyner AL, Costa RH, Bagchi S, Raychaudhuri P. FoxM1 in tumorigenicity of the neuroblastoma cells and renewal of the neural progenitors. Cancer Res. 2011;71:4292–302.View ArticleGoogle Scholar
- Xia LM, Huang WJ, Wang B, Liu M, Zhang Q, Yan W, Zhu Q, Luo M, Zhou ZZ, Tian DA. Transcriptional up-regulation of FoxM1 in response to hypoxia is mediated by HIF-1. J Cell Biochem. 2009;106:247–56.View ArticleGoogle Scholar
- Capasso M, Devoto M, Hou C, Asgharzadeh S, Glessner JT, Attiyeh EF, Mosse YP, Kim C, Diskin SJ, Cole KA, Bosse K, Diamond M, Laudenslager M, Winter C, Bradfield JP, Scott RH, Jagannathan J, Garris M, McConville C, London WB, Seeger RC, Grant SF, Li H, Rahman N, Rappaport E, Hakonarson H, Maris JM. Common variations in BARD1 influence susceptibility to high-risk neuroblastoma. Nat Genet. 2009;41:718–23.View ArticleGoogle Scholar
- Capasso M, Diskin SJ, Totaro F, Longo L, De Mariano M, Russo R, Cimmino F, Hakonarson H, Tonini GP, Devoto M, Maris JM, Iolascon A. Replication of GWAS-identified neuroblastoma risk loci strengthens the role of BARD1 and affirms the cumulative effect of genetic variations on disease susceptibility. Carcinogenesis. 2013;34:605–11.View ArticleGoogle Scholar
- Cimmino F, Avitabile M, Diskin SJ, Vaksman Z, Pignataro P, Formicola D, Cardinale A, Testori A, Koster J, de Torres C, Devoto M, Maris JM, Iolascon A, Capasso M. Fine mapping of 2q35 high-risk neuroblastoma locus reveals independent functional risk variants and suggests full-length BARD1 as tumor-suppressor. Int J Cancer. 2018;143(11):2828–37.View ArticleGoogle Scholar