Abstract
Francesco Maria Calabrese1,*, Rosanna Clima2,*, Piero Pignataro3,4, Vito Alessandro Lasorsa3,4, Michael D. Hogarty5, Aurora Castellano6, Massimo Conte7, Gian Paolo Tonini8, Achille Iolascon3,4, Giuseppe Gasparre2,#, Mario Capasso3,4,9,#
1Department of Biology, University of Bari “Aldo Moro”, Bari, Italy
2Department of Medical and Surgical Sciences-DIMEC, Medical Genetics Unit, University of Bologna, Bologna, Italy
3Department of Molecular Medicine and Medical Biotechnology, University of Naples Federico II, Naples, Italy
4CEINGE Biotecnologie Avanzate, Napoli, Italy
5Department of Pediatrics, Division of Oncology, Children’s Hospital of Philadelphia, Perelman School of Medicine at The University of Pennsylvania, Philadelphia, USA
6Paediatric Haematology/Oncology Department, IRCCS, Ospedale Pediatrico Bambino Gesù, Rome, Italy
7Oncology Unit, IRCCS Istituto Giannina Gaslini, Genova, Italy
8Pediatric Research Institute (IRP) – Fondazione Città della Speranza, Neuroblastoma Laboratory, Padua, Italy
9IRCCS SDN, Istituto di Ricerca Diagnostica e Nucleare, Naples, Italy
*These authors contributed equally to this work
#Co-last authors
Correspondence to:
Mario Capasso, email: [email protected]
Keywords: neuroblastoma, mitochondrial DNA mutations, WES, somatic mutations, germline mutations
Received: April 21, 2016 Accepted: June 09, 2016 Published: June 24, 2016
ABSTRACT
Background: Neuroblastoma, a tumor of the developing sympathetic nervous system, is a common childhood neoplasm that is often lethal. Mitochondrial DNA (mtDNA) mutations have been found in most tumors including neuroblastoma. We extracted mtDNA data from a cohort of neuroblastoma samples that had undergone Whole Exome Sequencing (WES) and also used snap-frozen samples in which mtDNA was entirely sequenced by Sanger technology. We next undertook the challenge of determining those mutations that are relevant to, or arisen during tumor development. The bioinformatics pipeline used to extract mitochondrial variants from matched tumor/blood samples was enriched by a set of filters inclusive of heteroplasmic fraction, nucleotide variability, and in silico prediction of pathogenicity.
Results: Our in silico multistep workflow applied both on WES and Sanger-sequenced neuroblastoma samples, allowed us to identify a limited burden of somatic and germline mitochondrial mutations with a potential pathogenic impact.
Conclusions: The few singleton germline and somatic mitochondrial mutations emerged, according to our in silico analysis, do not appear to impact on the development of neuroblastoma. Our findings are consistent with the hypothesis that most mitochondrial somatic mutations can be considered as ‘passengers’ and consequently have no discernible effect in this type of cancer.
BACKGROUND
Neuroblastoma (NB) originates in tissues of the sympathetic nervous system, typically within the adrenal medulla or paraspinal sympathetic ganglia. It is the most common neoplasm diagnosed in infancy and accounts for 10–15% of childhood cancer mortality [1]. Clinical and genetic risk factors are used to stratify patients into prognostic subgroups that differ in the predicted biological aggressiveness of disease [2]. Unfortunately, those children with severe clinical course and widespread metastases, often with unfavorable tumor genetics, are categorized as ‘high risk’ and have survival rates < 40% despite aggressive and intensive therapies [3]. The genetic basis of neuroblastoma has recently come into focus, as multiple genomic loci highly associated with sporadic neuroblastoma by genome-wide association studies (GWASs) have been discovered [4–11]. In addition, the tyrosine kinase receptor ALK, the most frequently mutated gene in familial neuroblastoma [12], somatically acquired amplification of MYCN [13] and hemizygous deletions of chromosomes 1p and 11q [4] are highly recurrent and associated with poor prognosis [10]. We and other groups have recently investigated the recurrence of nuclear somatic mutations in neuroblastoma by using diverse next generation approaches [14, 15]. Surprisingly, neuroblastomas have shown an exceeding low somatic nuclear mutation rate beyond those in ALK and ATRX loci but there are several less frequently mutated genes with significant implications in cancer pathways [15–17].
Moreover, very recently, a recurrent rearrangement that activates telomerase in high-risk neuroblastoma was reported. To date, the screening of neuroblastoma genomic mutations has mostly focused on the investigation of nuclear DNA. Only recently a paper has investigated the mitochondrial DNA (mtDNA) variants in primary and relapsed tumor from patients with neuroblastoma, concluding that nuclear and mitochondrial variants showed a concordant increase during tumor progression but the mutations with functional impact were rare [18].
Indeed, mutations can frequently occur in mtDNA and may contribute to tumorigenic processes by modifying metabolic, apoptotic or hypoxic mechanisms [19, 20]. The identification of clearly deleterious mtDNA mutations in cancer tissues - such as an intragenic deletion [21] or the common tRNALeu (UUR) A3243G MELAS mutation, or, even more so, very rare clearly pathogenic mutations in complex I [22] (for a review see [19, 23]) - underline the relevance of pathogenic mtDNA mutations in neoplastic transformation. In recent years a high number of point mutations, either germline or somatic, insertions and deletions in mtDNA have been identified in multiple cancers [24–28]. Moreover, a recent worldwide effort by the mitochondrial community has resulted in the creation of a consortium named MSeqDR [29], with the aim to collect, integrate, organize and critically analyze mitochondrial sequence data. Consequently, nuclear clinical databases together with genomic browsers may be used as a mold for the mitochondrial counterpart.
Hence, in order to provide an extensive characterization of the occurrence of mtDNA mutations in neuroblastoma, two different approaches were adopted. First, the MToolBox pipeline [30] was used to extract the mtDNA sequences from a training set of 26 whole exome sequencing (WES) samples. Secondly, mtDNA from 33 fresh-frozen neuroblastoma specimens was entirely sequenced by Sanger technology. Hence, for the first time, a detailed and stringent workflow was applied to filter the numerous variants found to depict the neuroblastoma mtDNA mutational burden, which may potentially impact on the disease course. By using a multi-parametric workflow for the prioritization of mitochondrial DNA variants of clinical interest, we finally selected a mutation pool conceivably pathogenic [31, 32].
RESULTS
Mitochondrial DNA extraction, variant detection and annotation
The mitochondrial variant characterization here presented was based on the application of a multistep workflow (Figure 1). Next Generation Sequencing (NGS) technologies provide a robust high-throughput batch of data to be checked for comprehensive mutation detection. This ought to also include an accurate and sensitive measurement of heteroplasmy. Because we consider unwise to define as potentially pathogenic or likely damaging most of the heteroplasmic variants emerging from high coverage read depth sequencing, we adopted a multi-parametric workflow previously developed [32] based on an exhaustive functional annotation through the mtDNA extraction pipeline MToolBox and including i) Macro Haplogroup Consensus Sequences to filter out fixed evolutionary variants, ii) rare or private variants prioritization, iii) the nucleotide variability as reported in HmtDB and iv) the disease score, based on several predictors of pathogenicity for non-synonymous variants [30].
Figure 1: Variant analysis and prioritization workflow: the main steps useful to the variant analysis and prioritization phases are based on the implementation of MToolBox pipeline with various filters and cutoffs. Disease score and nucleotide variability indexes are therefore accurately measured and evaluated. Finally, variants more likely contributing to defective phenotypes are selected.
Mitochondrial mean coverage read depth and mitochondrial assembled bases in WES datasets ranged from 6.25X to 2184.92X and from 19.72% to 100%, respectively (Supplementary_Table_1). A phylogeny quality check of samples was performed using MToolBox for haplogroup assignment. The quality score of the haplogroup predictions ranged between 85.7% and 100% (Supplementary_Table_1). When somatic/germline mitochondrial haplogroup assignment did not match, exception made for phylogenetic closely related branches (subgroups or higher-level connected hierarchy branches), paired samples were not considered.
Overall, as reported in Supplementary_Table_1, the clinical phenotype in our datasets included 40 High Risk (HR), 18 Low Risk (LR) and 1 Intermediate Risk (IR) patient, according to the Children’s Oncology Group (COG) classification system [2].
Because the mitochondrial assembled genomes and the depth of coverage for each samples were interconnected parameters useful in the quality sample estimate, in Supplementary Figure 1 the number of filtered mutations in both WES datasets was plotted against the relative mitochondrial mean of coverage depth in order to evaluate a potential correlation. Spearman’s rank correlation coefficient (rho and linked p-value) showed a non-linear correlation between the two considered variables in each dataset (see note in Supplementary Figure 1). The overall correlations were not significant, indicating that an increase in the mean coverage did not correspond with an increased filtered single nucleotide polymorphism (SNP) number.
Lack of association between specific mtDNA haplogroups and neuroblastoma
To investigate whether specific mtDNA haplogroups are associated to the occurrence of neuroblastoma, we compared the haplogroup frequencies observed in 43 Italian neuroblastoma samples within our cohort with those relative to 625 unrelated Italian healthy subjects, retrieved from the HmtDB public database [31] (Figure 2).
Figure 2: Macro-haplogroup frequency distribution in Italian neuroblastoma samples compared with the whole Italian healthy sample dataset retrieved from the HmtDB database shows the lack of association between specific mitochondrial haplogroups and disease occurrence.
A binary logistic regression was performed on the entire datasets and no significant p-value resulted, suggesting that no haplogroup might be significantly associated with the occurrence of this cancer.
WES germline variants
Pairwise comparison of blood and tumor mtDNA variants was performed in order to identify mtDNA germline variants. In 25/26 samples we identified a total of 362 variants, 244 of which were found against the three mitochondrial references used. Within this filtered variant batch, 142 (58.2%) defined haplogroups and were not considered in further pathogenic prioritization steps. Thirteen out of 102 were non-haplogroup defining events and showed non-overlapping confidence intervals (CI) of heteroplasmic fraction (HF) in tumor/blood pairs. Thus, we focused our attention on 89 variants with matching CI of HF (Supplementary_Table_2, sheet 1). The scatter plot shows the correlation between the HF values of germline variants, across the whole set of WES samples (Pearson’s product-moment correlation, r = 0.99, p-value = 2.2 × 10−16) (Figure 3A).
Figure 3: (A) Germline variant HF value correlation in blood and tumor. (B) The distribution of the heteroplasmic fractions (HFs) of germline variants shows how the majority of HFs (73%) had levels of heteroplasmy greater than 0.10.
Sixty-four of the remaining 89 variants exhibited an HF value greater than the fixed threshold (0.10) which ranged from 0.171 to 1, as reported in Figure 3B and Supplementary_Table_2 (sheet 2).
Twenty-three variants were synonymous (Supplementary_Table_2, sheet 3), 16 were non-synonymous (Supplementary_Table_2, sheet 4) and 25 belonged to the non-protein coding class (Supplementary_Table_2, sheet 5).
According to the disease score fixed threshold, higher than 0.4311, 7 out of 16 non-synonymous variants resulted to be potentially pathogenic in 6 samples out of 25 (24%) (Table 1).
Table 1: List of germiline mutations identified by WES on neuroblastoma samples
Sample | Clinical phenotype | Variant Allele | HF Blood/ Tumor | Locus | Nt Var | Aa Change | Disease Score | Mitomap Associated Disease(s) |
|---|---|---|---|---|---|---|---|---|
34C | HR | 11361C | 1/1 | MT-ND4 | 9.60 × 10−4 | M201T | 0.709 | |
1579N | LR | 14484C | 1/1 | MT-ND6 | 5.19 × 10−3 | M64V | 0.801 | LHON |
1C | HR | 9055A | 1/0.85 | MT-ATP6 | 2.34 × 10−1 | A177T | 0.533 | PD protective factor |
24C | HR | 11061G | 0.99/0.995 | MT-ND4 | 6.75 × 10−3 | S101C | 0.573 | |
8C | HR | 14180A | 0.962/1 | MT-ND6 | 2.18 × 10−2 | Y165F | 0.522 | |
8C | HR | 15665T | 0.615/1 | MT-CYB | 0 | L307F | 0.770 | |
NB08N | HR | 6419C | 0.171/0.082 | MT-CO1 | 0 | K172N | 0.782 |
Abbreviations: LR: Low Risk; HR: High Risk; HF: Heteroplasmic Fraction; Nt Var: nucleotide variability; LHON: Leber Hereditary Optic Neuropathy; PD: Parkinson Disease.
All reported non-synonymous mutations were detected as homoplasmic or nearly-homoplasmic in both blood and tumor, with the exception of the mutation m.6419A > C/MT-CO1 which exhibited a HF value close to the fixed threshold of 0.10. Although more than half of all samples harbored non-synonymous mutations in complex I gene subunits, the only one belonging to a low risk clinical phenotype (Italian sample 1579N) carried the deleterious homoplasmic mutation m.14484T > C (MT-ND6 Met64Val) also reported as associated to disease in Mitomap [33]. This mutation, strongly expected to predispose to a specific mitochondrial disease, was extensively published as one of the three most common mutations associated with Leber Hereditary Optic Neuropathy and responsible to induce a subtle reduction of complex I activity [19]. Moreover, the mutation m.9055G > A, found within the ATP6 gene in the sample 1C, was described to decrease Parkinson Disease risk [34]. All reported non-synonymous mutations had nucleotide variability ranging between 0 and 0.23. When we applied the nucleotide variability cutoff 0.0026 [32, 35], only three variants (m.11361T > C/MT-ND4, m.15665C > T/MT-CYB, m.6419A > C/MT-CO1) were classified as potentially pathogenic.
Among the 25 non protein-coding events, 20 mapped to the MT-DLOOP region of 15 individuals (Supplementary_Table_2, sheet 5). These batches included variants annotated as polymorphisms in Mitomap, as confirmed by their high nucleotide variability score and the reported allele frequency in the 1000 Genomes database. All 20 MT-DLOOP variants were homoplasmic (or nearly homoplasmic) and all defined a different haplogroup than the assigned one. Among non-protein coding variants, two resided within the rRNA12S locus of the same sample and three in tRNA loci. The secondary structure of mutated rRNA12S showed no modification and a nearly unvaried minimum free energy (data not shown). Three variants belonged to tRNA loci, namely the m.14687A > G/MT-TE, the m.15946C > T/MT-TT and the m.10463T > C/MT-TR. PhyloP and PhastCons tRNA conservation scores, consulted to check variant position preservation during evolution and functional constraints, were equal to 1.56 and 0.82, −0.98 and 0, −1.55 and 0 in glutamic acid, threonine and arginine tRNAs, respectively. These different values pointed out how the position m.14687 of tRNA glutamic acid resulted to be phylogenetically more conserved than the other two. All germline tRNA variants contributed to define other haplogroups, somewhat diminishing the possibility of a pathogenetic role.
Briefly, although the germline defined variants were overall 362, no one of these appeared to be potentially pathogenic within the non-protein-coding subset; only 7 (2%) complied with all prioritization criteria, mapped in protein-coding genes, and were therefore selected as potentially pathogenic mutations (Table 1).
WES somatic variants
Those events uniquely detected in tumor samples were annotated as somatic variants. To assess this status, a further check was performed on each predicted germline paired variant, looking for those sites not detectable due to a coverage lower than the minimum of 5 in MToolBox. We identified 119 false positive somatic mutations (Supplementary_Table_3, sheet 1). The analysis was thus centered on a total of 60 somatic variants, all recognized with respect to the three mitochondrial references used (Supplementary_Table_3, sheet 2). Of these, 59 variants did not contribute to define the sample haplogroup (Supplementary_Table_3, sheet 3). The majority of them (89.8%) showed a HF lower than 0.10 and were consequently discarded (Figure 4). The remaining 6 variants had low HF (0.104 – 0.211) and the nucleotide variability did not exceed the value of 0.011 (Supplementary_Table_3, sheet 4); one was synonymous (16.7%) (m.8239C > T/MT-CO2), and five belonged to the non-synonymous variant class (Supplementary_Table_3, sheet 5).
Figure 4: Distribution of the heteroplasmic fractions (HFs) of somatic variants.
Among the non-synonymous variants, 4 were recognized as potentially pathogenic (Table 2) in as many different samples. The mutation m.11090A > C/MT-ND4 with a disease score of 0.85 co-occurred in two samples (2221T and 1940T). The 34S sample had a mutation in the MT-CO3 gene with the highest HF value (0.208) and a relative null nucleotide variability value (Table 2). All potentially deleterious WES somatic mutations did not contribute to define other haplogroups and were not annotated in Mitomap or OMIM in association to disease.
Table 2: List of somatic mutations identified by WES on neuroblastoma samples
Sample | Clinical phenotype | Variant Allele | HF | Locus | Nt Var | Aa change | Disease Score |
|---|---|---|---|---|---|---|---|
34S | HR | 9625T | 0.208 | MT-CO3 | 0 | S140L | 0.857 |
NB11PT | HR | 10306C | 0.125 | MT-ND3 | 0 | N83T | 0.503 |
2221T | HR | 11090C | 0.122 | MT-ND4 | 5.10 × 10−4 | T111P | 0.850 |
1940T | HR | 11090C | 0.104 | MT-ND4 | 5.10 × 10−4 | T111P | 0.850 |
Abbreviations: HR: High Risk; HF: Heteroplasmic Fraction; Nt Var: nucleotide variability.
Overall, all non-synonymous variants showed nucleotide variability and disease score respectively lower and higher than the fixed thresholds.
Mitochondrial DNA mutations in snap-frozen neuroblastoma samples
Full mtDNA sequencing was also performed by Sanger technology in 33 snap-frozen tumor biopsies (Supplementary_Table_1). The corresponding blood-extracted DNA allowed us to estimate the germline or somatic nature of potentially pathogenic mutations in our set. A total of 476 mtDNA events were identified (Supplementary_Table_4, sheet 1). Of this first batch, 287 were recognized by all the three mitochondrial references (Supplementary_Table_4, sheet 2). The variants contributing to haplogroup assignment were 191 and were discarded from the subsequent analyses. The analysis was centered on 96 non-haplogroup defining variants (Supplementary_Table_4, sheet 3), 78 of which mapped within protein-coding regions (Supplementary_Table_4, sheet 4) and 18 in rRNAs and MT-DLOOP loci (Supplementary_Table_4, sheet 5). The non-protein-coding variants displayed variability values equal to or lower than 0.031. Specifically, 2 resided within the regulatory MT-DLOOP region and were not previously reported in Mitomap. Moreover, 7 patients harbored variants in tRNA genes. Four did not define other haplogroups and had the lowest nucleotide variability score (Table 3). When we looked at the tRNA secondary structure, we noticed that the somatic m.7506G > A variant in the serine tRNA (MT-TS1) locus mapped in the D arm stem and was previously reported as associated with progressive external ophthalmoplegia and hearing loss, in heteroplasmy [36] (Table 3). Two germline variants mapped on tRNA-Thr (MT-TT) both in the stem of the T arm: the first m.15943T > C did not impact on secondary structure, the second one (m.15944T > C) is annotated in MAMIT-tRNA [37] as a polymorphism (Table 3). We found the germline m.10023C > G/MT-TG mutation, whose position is in a highly evolutionarily conserved site, i.e. the anticodon stem-loop (Table 3). The m.10023C > G mutation would cause the carrying of the amino acid arginine instead of glycine. The two amino acids belong to different chemical classes; the non-polar charged glycine would therefore be substituted with a hydrophilic polar charged one. We therefore deemed this one mutation to be the only one potentially interesting from the clinical point of view, although the lack of data on the segregation of the same variant in the individual’s mother does not allow us to conclude on a potential predisposing effect on neuroblastoma of this mutation.
Table 3: List of mutations identified in tRNAs by Sanger sequencing on snap-frozen neuroblastoma samples
Sample | Clinical phenotype | Variant Allele | Status | HF1 | Locus | Nt Var | Mitomap Associated Disease(s) |
|---|---|---|---|---|---|---|---|
1360 | HR | 10023G | germline | MT-TG | 0 | ||
1506 | LR | 7506A | somatic | MT-TS1 | 0 | PEO with hearing loss | |
1493 | HR | 15943C | germline | MT-TT | 1.0 × 10−3 | ||
1360 | HR | 15944C | germline | MT-TT | 2.0 × 10−3 |
Abbreviations; HR: High Risk; LR: Low Risk; Nt Var: nucleotide variability; PEO: Progressive External Ophthalmoplegia. 1Heteroplasmic fraction is not available for mutations identified with Sanger sequencing. Here only heteroplasmic (+) vs. homoplasmic (−) status is reported.
Ribosomal RNA genes were also found to harbor 7 substitutions. Four of them were in the MT-RNR2 and 3 in the MT-RNR1 gene, in seven different samples (Supplementary_Table_4, sheet 5). Prediction of the secondary structure of the rRNA mutated molecules showed no evident structural changes for all substitutions compared to the wild-type conformation (data not shown).
Twenty-six non-synonymous variants (Supplementary_Table_4; sheet 6) were detected in 17 out of 33 samples (51.5%), 10 of which were predicted to impair protein function by a disease score higher than 0.4311. This number was further reduced to 6 potentially pathogenic variants, when the nucleotide variability cutoff filter was applied (Supplementary_Table_4, sheet 7). All of these considered variants were exclusive to each sample. In addition, one frame-shift event was also detected (Supplementary_Table_4, sheet 7).
Upon assessing the germline/somatic status of 7 of these potential pathogenic mutations in protein coding regions, we found 2 germline and 5 somatic mutations (Table 4).
Table 4: List of prioritized mutations identified by Sanger sequencing on snap-frozen neuroblastoma samples
Sample | Clinical phenotype | Variant Allele | Status | HF1 | Locus | Nt Var | AA Change | Disease Score |
|---|---|---|---|---|---|---|---|---|
1332 | HR | 15284G | germline | − | MT-CYB | 0 | T180A | 0.626 |
1360 | HR | 12014T | germiline | + | MT-ND4 | 2.0 × 10−3 | L419F | 0.524 |
1090 | LR | 7402A | somatic | − | MT-CO1 | 0 | P500Q | 0.748 |
1506 | LR | 8441A | somatic | − | MT-ATP8 | 1.0 × 10−3 | L26M | 0.784 |
9056A | somatic | − | MT-ATP6 | 1.0 × 10−3 | A177D | 0.729 | ||
1047 | HR | 14945A | somatic | − | MT-CYB | 1.0 × 10−3 | A67T | 0.678 |
1144 | HR | 11915d | somatic | − | MT-ND4 | 0 | frameshift |
Abbreviations: HR: High Risk; HF: Heteroplasmic Fraction; Nt Var: nucleotide variability. 1Heteroplasmic fraction is not available for mutations identified with Sanger sequencing. Here only heteroplasmic (+) vs. homoplasmic (−) status is reported.
None of the germline and somatic mutations were found to be associated to disease in Mitomap or OMIM. Sample 1360 harbored the only one mutation found in heteroplasmy (C/T) in the blood (m.12014C > T/MT-ND4). The remaining 2 germline mutations mapped in the gene MT-ND4 encoding for a complex I subunit (Table 4).
All 5 reported somatic mutations had a high disease score ranging between 0.678 and 0.784 (Table 4). One frameshift mutation occurring in the MT-ND4 gene was detected (Table 4).
DISCUSSION
In this study, we intended to ascertain whether mitochondrial DNA variants might be linked to neuroblastoma. The MToolBox pipeline and a multi-parametric workflow, both recently developed [30, 32], were used to investigate and prioritize the germline and somatic mutations in two WES datasets of tumors and blood pairs of 26 individuals. Moreover, full mtDNA from 33 snap-frozen neuroblastoma samples was sequenced with Sanger method and analyzed with MToolBox as well. A similar strategy was previously exploited in an analogous fashion by our group on glioblastoma multiforme [38]. The identification and the understanding of genetic events that drive cancer is a crucial phase to improve diagnostic, prognostic and therapeutic strategies. Hereby, we here intended to explore the mitochondrial variation and its possible involvement in neuroblastoma.
Neuroblastoma, like most other solid tumors, shows a high glucose uptake [39, 40], together with a high rate of lactic acid production and a low rate of oxygen consumption, reflecting the switch from mitochondrial OXPHOS to glycolysis [41]. Moreover, a down-regulation of all components of the aerobic mitochondrial energy metabolism without affecting mitochondrial mass was recently observed [42]. Recently, an in vitro study on the neuroblastoma cell lines LA-N-1, IMR-32, LS and SK-N-SH, showing an increased oxidative stress, a reduced lactate dehydrogenase (LDH) enzyme activity and reduced lactate production after the use of the antiprotozoal drug nifurtimox, supports the previous finding [43]. Moreover, it is becoming clearer how many oncogenic key signaling pathways converge to adapt tumor cell metabolism in order to support cancer growth and survival. Indeed, the glycolytic shift in the presence of oxygen (Warburg effect) has been inferred in many cancers, including neuroblastoma [39, 41, 44]. In the last decade several groups have investigated changes of the OXPHOS system in cancer. It has been proven that, depending on the tumor type, cancer cells are frequently characterized by either a generally low amount of mitochondria and OXPHOS [45], a specific defect of one OXPHOS complex [46–49] or a combined low level of OXPHOS complexes [50]. There are strong evidences that the Warburg effect in neuroblastoma is not caused by a low mitochondrial mass per cell [51]. Thus, defects of mitochondrial functions may be relevant modifiers events in neuroblastoma, from which derived our need to investigate whether this cancer type may represent a paradigm in which the Warburg effect is due to mtDNA mutations. Nonetheless, our data do not point to a major involvement of mtDNA mutations in neuroblastoma genesis or progression, as the burden of the genetic lesions that pass stringent prioritization filters for a potential pathogenic role was limited.
Overall, we found very few somatic mutations, in agreement with recent high-throughput screenings of nuclear DNAs which have demonstrated the relative paucity of recurrent somatic mutations in neuroblastoma and highlighted the difficulty to develop therapeutic strategies, relying on frequently altered oncogenic drivers [14, 15].
Specifically for mtDNA, one potential explanation for this is that these tumors arise in very young individuals, whereas accumulation of mtDNA mutations is a phenomenon typical of aging tissues. It is also true, however, that it is becoming increasingly evident that most mtDNA mutations, particularly if detected in low heteroplasmies, may result from intrinsic mtDNA replication mechanisms and are to be considered passengers events during tumorigenesis, few of which only becoming fixed in the quickly-evolving cancer cell population [52]. With respect to the few potentially pathogenic mutations found, it is worth to mention that, although aware of the low statistical power to attempt a correlation, we found no association of the former with MYCN ploidy, based on the rationale that the proto-oncogene MYCN is indirectly involved in regulating the switch off of mitochondrial function via a positive regulation of aerobic glycolysis and conversion of pyruvate to lactate [43, 53]. Its amplification, therefore, besides being associated to a poorer prognosis in neuroblastoma, may favor a relaxed selection of mtDNA variants. This does not appear to be the case, based on the limited data we gathered in our cohorts, albeit further studies are warranted on this issue.
Overall, apart from the environmental factors acting over mtDNA mutations which may easily concur to this multifactorial neoplasm, our in silico data, when merged with the other evidences present in the scientific literature on the metabolic features of neuroblastoma, suggest a mechanism of OXPHOS reduction likely independent from the occurrence of mtDNA mutations.
MATERIALS AND METHODS
Ethics statement
The present study was formally conducted on two different datasets, one composed of Italian and one of American neuroblastoma patients. Data relative to the first one, obtained both with Sanger and next generation sequencing (NGS) technology, were approved by the Ethics Committee of the Medical University of Naples; the written informed consent was obtained by all children’s legal guardians. The other batch of analyzed WES, stored in the European Genome-phenome Archive (EGA) database, was used with the approval obtained from The Children’s Hospital of Philadelphia Institutional Review Board. This second dataset is part of an enlarged neuroblastoma study published by Sausen in 2013 [54].
WES Children’s Oncology Group (COG) samples
Patient details and any associated reference for the 16 matched primary tumor/normal samples belonging to the COG clinical trials cooperative group biobank here used, are retrievable in the online version of the previously cited paper (REF) and resulted as part of the EGA project entry registered under Acc. No EGAS00001000369.
WES and snap-frozen DNA neuroblastoma samples
In-house WES data of neuroblastoma samples were used to retrieve the mtDNA sequences. DNA samples for Sanger sequencing were obtained from the BioBank at IRCCS AOU San Martino-IST. Agilent SureSelect Target Enrichment System 50 Mb (Agilent Technologies, Santa Clara, CA) was used for exome capture (10 normal/tumor DNA) according to the manufacturers’ protocol. Captured DNAs were subjected to massively parallel sequencing using Illumina HiSeq 2000 (Illumina Inc., San Diego, CA) with 100 bp paired-end reads, while 33 DNA samples were sequenced by using Sanger sequencing technology.
mtDNA amplification and sequencing (sanger method)
Single PCR volumes were optimized in order to allow a better performance amplification on mtDNA and were set as follows: 25 μl in order to allow a better performance amplification (HiFi HotStart DNA Polymerase - KAPABIOSYSTEMS), 3.5 μl primer forward + reverse each one concentrated 10 μM, 6 μl DNA (30 ng tot.) and 3 μl H2O. 46 overlapping amplicons in order to cover full-length mtDNA (available upon request) were used. Cycle amplification were made of 3′ to 95°C, 15″ to 95°C, 5″ to 60°C and 5″ to 72°C (× 40 times), finally 10′ to 72°C. The amplified amplicons were purified and sequenced by using a 3730 DNA Analyzer (Applied Biosystem).
Haplogroup frequency statistics
A set of 625 Italian healthy samples retrieved from the human mitochondrial public database HmtDB [31] was used as haplogroup control group and was tested against our Italian neuroblastoma dataset.
Statistical analyses were performed using R environment, version 3.1.2, and statistical significance was established at p-value ≤ 0.05. Binary logistic-regression analysis was used. In order to test a potential haplogroup contribution to neuroblastoma etiology or predisposition, a logistic regression of all categories at once was performed. We analyzed the haplogroup categories in patients and controls using the common haplogroup H as the reference level. Spearman rank test was used to assess the correlation between mitochondrial mean coverage read depth and the number of filtered variants. Pearson product-moment correlation coefficient, a measure of the strength and direction of association that exists between two continuous variables, was applied.
Retrieving of mtDNA reads and prioritization criteria
All WES dataset fastq and fasta Sanger sequence files were analyzed by using the MToolBox pipeline [30], which helps filtering and prioritizing pathogenic variants.
For each tumor/blood samples, three parameter discrepancies resulting from MtoolBox, i.e. mitochondrial read depth, assembled bases and haplogroup assignment were considered; if these values for a tumor/normal pair were not comparable, samples were discarded from the analysis. Moreover, because the read depth threshold required by MToolBox filtered out variants with a depth lower than 5×, the germline counterpart in the mtDNAassembly-table output was checked, in order to manually assess, by inspecting the coverage, the somatic status for each detected variant. Further, we considered as potentially pathogenic variants only those recognized against the mitochondrial reference sequences (revised Cambridge Reference Sequence (rCRS), Recostructed Sapiens Reference Sequence (RSRS) and Macro Haplogroup Consensus Sequence (MHCS)) which occurred in non haplogroup-defining sites, which featured a heteroplasmy fraction greater than 0.10, featuring nucleotide variability lower than nucleotide variability cutoff (0.0026) and having a disease score above the disease score threshold (0.4311). We also checked the Heteroplasmy Fraction (HF) confidence interval related concordance in variants shared between tumor/normal-paired samples. For the whole set of selected variants the nucleotide variability value, both for somatic and germline samples, was estimated [35]. Intuitively, mitochondrial low variability values, directly linked to positions not prone to vary, were associated more likely to potentially pathogenic mutations [55]. To evaluate if the number of mitochondrial variant positions (found versus all the three reference sequences and singularly – without any other applied filter) and MYCN amplification status correlated, a Wilcoxon test (R environment) was performed.
The evaluation of tRNA variants was based on the phastCons score [56], which takes into account neighboring bases and reports the probability that each nucleotide belongs to a conserved element, and on the phyloP algorithm [57], useful in evaluating signatures of selection at particular nucleotides, ignoring neighboring bases in its calculation and producing a p-value. The second algorithm measures both acceleration (faster evolution than expected under neutral drift) and conservation (slower than expected evolution).
The prediction of RNA secondary structures was performed by RNA Folding Form (version 2.3 energies) software, using the set of default parameters [58].
Heteroplasmy fraction threshold
Because of the uncertain potential degree of contamination affecting the tumor tissue, which is not measurable by analyzing the non-tumor counterpart, a heteroplasmy underestimation ought here to be considered. Further, another limitation that may affect the mitochondrial coverage and therefore the prioritization process, deals with the enrichment technology used. As noticed in human exome projects, different enrichment technologies showed different off-target read distributions [59]. The highest mitochondrial read coverage distributions were observed in association with lower density platforms. Both the WES sequencing used in this paper are based on SureSelect Agilent enrichment technology which showed an intermediately mitochondrial coverage read depth distribution together with a corresponding reconstructed genome percentage in terms of reconstructed contigs. The SuresSelect ranks between other two powerful used platforms: the Illumina TrueSeq and Nimblegen SeqCap EZ-exome platform. Thus, after considering both the uncertain level of tumor tissue contamination and the off-target mitochondrial spectrum platform dependence, we forthwith set the HF threshold to 10%. A similar HF threshold was used for other heterogeneous tumor like glioblastoma [38].
CONCLUSIONS
Although a high number of mitochondrial variants was detected, the large majority of them did not pass the prioritization filters, as the WES detecting capability can be considered a double edge sword, which sums up a huge number of mutations which we opportunely filtered in order to recognize those of potential clinical interest. Therefore, a relatively limited burden of pathogenic mutations is indeed carried by neuroblastoma suggesting that they may be passenger events rather than potent modifiers in cancer progression in this cancer type. More efforts are warranted to understand if mtDNA contributes at all to neuroblastoma aggressive behavior, its progression and its metabolic plasticity.
ACKNOWLEDGMENTS AND FUNDING
The work was supported in part by Associazione Italiana per la Ricerca sul Cancro – AIRC grant n.10537 to MC, and grant IG14242 to GG; by Ministero della Salute (GR-2011-02348722) to MC and GG; Fondazione Italiana per la Lotta al Neuroblastoma (MC, GPT and LL) and Associazione Oncologia Pediatrica e Neuroblastoma (MC). We thank INFN at Physic department (University of Bari) for computational computing resources availability.
CONFLICTS OF INTEREST
The authors declare no conflicts of interest.
Authors’ contributions
FMC, RC and GG analyzed data and wrote the manuscript. PP performed the wet bench amplification and sequencing. AL took care of data downloads and transfer protocols. MDH gathered and provided WES American samples BAM files together with the related clinical data. AC, MC, GPT provided NBL Italian patients clinical data and follow the project development. AI, GG and MC coordinated and supervised the whole project. All authors read and approved the final manuscript.
REFERENCES
1. Brodeur GM. Neuroblastoma: biological insights into a clinical enigma. Nat Rev Cancer. 2003; 3:203–16. doi: 10.1038/nrc1014.
2. Cohn SL, Pearson ADJ, London WB, Monclair T, Ambros PF, Brodeur GM, Faldum A, Hero B, Iehara T, Machin D, Mosseri V, Simon T, Garaventa A, et al. The International Neuroblastoma Risk Group (INRG) classification system: an INRG Task Force report. J Clin Oncol Off J Am Soc Clin Oncol. 2009; 27:289–97. doi: 10.1200/JCO.2008.16.6785.
3. Maris JM. Recent advances in neuroblastoma. N Engl J Med. 2010; 362:2202–11. doi: 10.1056/NEJMra0804577.
4. Maris JM, Mosse YP, Bradfield JP, Hou C, Monni S, Scott RH, Asgharzadeh S, Attiyeh EF, Diskin SJ, Laudenslager M, Winter C, Cole KA, Glessner JT, et al. Chromosome 6p22 locus associated with clinically aggressive neuroblastoma. N Engl J Med. 2008; 358:2585–93. doi: 10.1056/NEJMoa0708698.
5. 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, et al. Common variations in BARD1 influence susceptibility to high-risk neuroblastoma. Nat Genet. 2009; 41:718–23. doi: 10.1038/ng.374.
6. Nguyen LB, Diskin SJ, Capasso M, Wang K, Diamond MA, Glessner J, Kim C, Attiyeh EF, Mosse YP, Cole K, Iolascon A, Devoto M, Hakonarson H, et al. Phenotype restricted genome-wide association study using a gene-centric approach identifies three low-risk neuroblastoma susceptibility Loci. PLoS Genet. 2011; 7:e1002026. doi: 10.1371/journal.pgen.1002026.
7. Wang K, Diskin SJ, Zhang H, Attiyeh EF, Winter C, Hou C, Schnepp RW, Diamond M, Bosse K, Mayes PA, Glessner J, Kim C, Frackelton E, et al. Integrative genomics identifies LMO1 as a neuroblastoma oncogene. Nature. 2011; 469:216–20. doi: 10.1038/nature09609.
8. Diskin SJ, Capasso M, Schnepp RW, Cole KA, Attiyeh EF, Hou C, Diamond M, Carpenter EL, Winter C, Lee H, Jagannathan J, Latorre V, Iolascon A, et al. Common variation at 6q16 within HACE1 and LIN28B influences susceptibility to neuroblastoma. Nat Genet. 2012; 44:1126–30. doi: 10.1038/ng.2387.
9. Diskin SJ, Capasso M, Diamond M, Oldridge DA, Conkrite K, Bosse KR, Russell MR, Iolascon A, Hakonarson H, Devoto M, Maris JM. Rare variants in TP53 and susceptibility to neuroblastoma. J Natl Cancer Inst. 2014; 106:dju047. doi: 10.1093/jnci/dju047.
10. 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. doi: 10.1093/carcin/bgs380.
11. Capasso M, Diskin S, Cimmino F, Acierno G, Totaro F, Petrosino G, Pezone L, Diamond M, McDaniel L, Hakonarson H, Iolascon A, Devoto M, Maris JM. Common genetic variants in NEFL influence gene expression and neuroblastoma risk. Cancer Res. 2014; 74:6913–24. doi: 10.1158/0008-5472.CAN-14-0431.
12. Mossé YP, Laudenslager M, Longo L, Cole KA, Wood A, Attiyeh EF, Laquaglia MJ, Sennett R, Lynch JE, Perri P, Laureys G, Speleman F, Kim C, et al. Identification of ALK as a major familial neuroblastoma predisposition gene. Nature. 2008; 455:930–5. doi: 10.1038/nature07261.
13. Schwab M. Amplified MYCN in human neuroblastoma: paradigm for the translation of molecular genetics to clinical oncology. Ann N Y Acad Sci. 2002; 963:63–73.
14. Schulte JH, Eggert A. Neuroblastoma. Crit Rev Oncog. 2015; 20:245–70.
15. 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, et al. Exome and deep sequencing of clinically aggressive neuroblastoma reveal somatic mutations that affect key pathways involved in cancer progression. Oncotarget. 2016; . doi: 10.18632/oncotarget.8187.
16. Eleveld TF, Oldridge DA, Bernard V, Koster J, Daage LC, Diskin SJ, Schild L, Bentahar NB, Bellini A, Chicard M, Lapouble E, Combaret V, Legoix-Né P, et al. Relapsed neuroblastomas show frequent RAS-MAPK pathway mutations. Nat Genet. 2015; 47:864–71. doi: 10.1038/ng.3333.
17. Schramm A, Köster J, Assenov Y, Althoff K, Peifer M, Mahlow E, Odersky A, Beisser D, Ernst C, Henssen AG, Stephan H, Schröder C, Heukamp L, et al. Mutational dynamics between primary and relapse neuroblastomas. Nat Genet. 2015; 47: 872–7. doi: 10.1038/ng.3349.
18. Riehl LM, Schulte JH, Mulaw MA, Dahlhaus M, Fischer M, Schramm A, Eggert A, Debatin K-M, Beltinger C. The mitochondrial genetic landscape in neuroblastoma from tumor initiation to relapse. Oncotarget. 2016; 7:6620–5. doi: 10.18632/oncotarget.6776.
19. Iommarini L, Calvaruso MA, Kurelac I, Gasparre G, Porcelli AM. Complex I impairment in mitochondrial diseases and cancer: parallel roads leading to different outcomes. Int J Biochem Cell Biol. 2013; 45:47–63. doi: 10.1016/j.biocel.2012.05.016.
20. Wallace DC. Mitochondria and cancer. Nat Rev Cancer. 2012; 12: 685–98. doi: 10.1038/nrc3365.
21. Horton TM, Petros JA, Heddi A, Shoffner J, Kaufman AE, Graham SD, Gramlich T, Wallace DC. Novel mitochondrial DNA deletion found in a renal cell carcinoma. Genes Chromosomes Cancer. 1996; 15:95–101. doi: 10.1002/(SICI)1098-2264(199602)15:2<95::AID-GCC3>3.0.CO;2-Z.
22. Meierhofer D, Mayr JA, Fink K, Schmeller N, Kofler B, Sperl W. Mitochondrial DNA mutations in renal cell carcinomas revealed no general impact on energy metabolism. Br J Cancer. 2006; 94:268–74. doi: 10.1038/sj.bjc.6602929.
23. Gasparre G, Porcelli AM, Lenaz G, Romeo G. Relevance of mitochondrial genetics and metabolism in cancer development. Cold Spring Harb Perspect Biol. 2013; 5. doi: 10.1101/cshperspect.a011411.
24. Bai R-K, Leal SM, Covarrubias D, Liu A, Wong L-JC. Mitochondrial genetic background modifies breast cancer risk. Cancer Res. 2007; 67:4687–94. doi: 10.1158/0008-5472.CAN-06-3554.
25. Booker LM, Habermacher GM, Jessie BC, Sun QC, Baumann AK, Amin M, Lim SD, Fernandez-Golarz C, Lyles RH, Brown MD, Marshall FF, Petros JA. North American white mitochondrial haplogroups in prostate and renal cancer. J Urol. 2006; 175:468–472–473. doi: 10.1016/S0022-5347(05)00163-1.
26. Gómez-Zaera M, Abril J, González L, Aguiló F, Condom E, Nadal M, Nunes V. Identification of somatic and germline mitochondrial DNA sequence variants in prostate cancer patients. Mutat Res. 2006; 595:42–51. doi: 10.1016/j. mrfmmm.2005.10.012.
27. Kloss-Brandstätter A, Schäfer G, Erhart G, Hüttenhofer A, Coassin S, Seifarth C, Summerer M, Bektic J, Klocker H, Kronenberg F. Somatic mutations throughout the entire mitochondrial genome are associated with elevated PSA levels in prostate cancer patients. Am J Hum Genet. 2010; 87: 802–12. doi: 10.1016/j.ajhg.2010.11.001.
28. Lam ET, Bracci PM, Holly EA, Chu C, Poon A, Wan E, White K, Kwok P-Y, Pawlikowska L, Tranah GJ. Mitochondrial DNA sequence variation and risk of pancreatic cancer. Cancer Res. 2012; 72: 686–95. doi: 10.1158/0008-5472.CAN-11-1682.
29. Falk MJ, Shen L, Gonzalez M, Leipzig J, Lott MT, Stassen APM, Diroma MA, Navarro-Gomez D, Yeske P, Bai R, Boles RG, Brilhante V, Ralph D, et al. Mitochondrial Disease Sequence Data Resource (MSeqDR): a global grass-roots consortium to facilitate deposition, curation, annotation, and integrated analysis of genomic data for the mitochondrial disease clinical and research communities. Mol Genet Metab. 2015; 114: 388–96. doi: 10.1016/j.ymgme.2014.11.016.
30. Calabrese C, Simone D, Diroma MA, Santorsola M, Guttà C, Gasparre G, Picardi E, Pesole G, Attimonelli M. MToolBox: a highly automated pipeline for heteroplasmy annotation and prioritization analysis of human mitochondrial variants in high-throughput sequencing. Bioinformatics. 2014; 30: 3115–7. doi: 10.1093/bioinformatics/btu483.
31. Rubino F, Piredda R, Calabrese FM, Simone D, Lang M, Calabrese C, Petruzzella V, Tommaseo-Ponzetta M, Gasparre G, Attimonelli M. HmtDB, a genomic resource for mitochondrion-based human variability studies. Nucleic Acids Res. 2012; 40:D1150–9. doi: 10.1093/nar/gkr1086.
32. Santorsola M, Calabrese C, Girolimetti G, Diroma MA, Gasparre G, Attimonelli M. A multi-parametric workflow for the prioritization of mitochondrial DNA variants of clinical interest. Hum Genet. 2016; 135:121–36. doi: 10.1007/s00439-015-1615-9.
33. Lott MT, Derbeneva O, Chalkia D, Sarmady M, Procaccio V, Wallace DC. mtDNA variation and analysis using MITOMAP and MITOMASTER. [cited 2016 May 17]. Available 2016 May 17, from http://www.mitomap.org.
34. van der Walt JM, Nicodemus KK, Martin ER, Scott WK, Nance MA, Watts RL, Hubble JP, Haines JL, Koller WC, Lyons K, Pahwa R, Stern MB, Colcher A, et al. Mitochondrial polymorphisms significantly reduce the risk of Parkinson disease. Am J Hum Genet. 2003; 72:804–11.
35. Pesole G, Saccone C. A novel method for estimating substitution rate variation among sites in a large dataset of homologous DNA sequences. Genetics. 2001; 157:859–65.
36. Cardaioli E, Da Pozzo P, Gallus GN, Malandrini A, Gambelli S, Gaudiano C, Malfatti E, Viscomi C, Zicari E, Berti G, Serni G, Dotti MT, Federico A. A novel heteroplasmic tRNA(Ser(UCN)) mtDNA point mutation associated with progressive external ophthalmoplegia and hearing loss. Neuromuscul Disord NMD. 2007; 17:681–3. doi: 10.1016/j. nmd.2007.05.001.
37. Pütz J, Dupuis B, Sissler M, Florentz C. Mamit-tRNA, a database of mammalian mitochondrial tRNA primary and secondary structures. RNA. 2007; 13:1184–90. doi: 10.1261/rna.588407.
38. Vidone M, Clima R, Santorsola M, Calabrese C, Girolimetti G, Kurelac I, Amato LB, Iommarini L, Trevisan E, Leone M, Soffietti R, Morra I, Faccani G, et al. A comprehensive characterization of mitochondrial DNA mutations in glioblastoma multiforme. Int J Biochem Cell Biol. 2015; 63:46–54. doi: 10.1016/j.biocel.2015.01.027.
39. Kushner BH, Yeung HW, Larson SM, Kramer K, Cheung NK. Extending positron emission tomography scan utility to high-risk neuroblastoma: fluorine-18 fluorodeoxyglucose positron emission tomography as sole imaging modality in follow-up of patients. J Clin Oncol Off J Am Soc Clin Oncol. 2001; 19:3397–405.
40. Freebody J, Wegner EA, Rossleigh MA. 2-deoxy-2-((18)F) fluoro-D-glucose positron emission tomography/computed tomography imaging in paediatric oncology. World J Radiol. 2014; 6:741–55. doi: 10.4329/wjr.v6.i10.741.
41. Aminzadeh S, Vidali S, Sperl W, Kofler B, Feichtinger RG. Energy metabolism in neuroblastoma and Wilms tumor. Transl Pediatr. 2015; 4:20–32. doi: 10.3978/j.issn.2224-4336.2015.01.04.
42. Feichtinger RG, Zimmermann F, Mayr JA, Neureiter D, Hauser-Kronberger C, Schilling FH, Jones N, Sperl W, Kofler B. Low aerobic mitochondrial energy metabolism in poorly- or undifferentiated neuroblastoma. BMC Cancer. 2010; 10: 149. doi: 10.1186/1471-2407-10-149.
43. Cabanillas Stanchi KM, Bruchelt G, Handgretinger R, Holzer U. Nifurtimox reduces N-Myc expression and aerobic glycolysis in neuroblastoma. Cancer Biol Ther. 2015; 16:1353–63. doi: 10.1080/15384047.2015.1070987.
44. Smith DJ, Cossins LR, Hatzinisiriou I, Haber M, Nagley P. Lack of correlation between MYCN expression and the Warburg effect in neuroblastoma cell lines. BMC Cancer. 2008; 8:259. doi: 10.1186/1471-2407-8-259.
45. Meierhofer D, Mayr JA, Foetschl U, Berger A, Fink K, Schmeller N, Hacker GW, Hauser-Kronberger C, Kofler B, Sperl W. Decrease of mitochondrial DNA content and energy metabolism in renal cell carcinoma. Carcinogenesis. 2004; 25:1005–10. doi: 10.1093/carcin/bgh104.
46. Zimmermann FA, Mayr JA, Feichtinger R, Neureiter D, Lechner R, Koegler C, Ratschek M, Rusmir H, Sargsyan K, Sperl W, Kofler B. Respiratory chain complex I is a mitochondrial tumor suppressor of oncocytic tumors. Front Biosci Elite Ed. 2011; 3:315–25.
47. Calabrese C, Iommarini L, Kurelac I, Calvaruso MA, Capristo M, Lollini P-L, Nanni P, Bergamini C, Nicoletti G, Giovanni CD, Ghelli A, Giorgio V, Caratozzolo MF, et al. Respiratory complex I is essential to induce a Warburg profile in mitochondria-defective tumor cells. Cancer Metab. 2013; 1:11. doi: 10.1186/2049-3002-1-11.
48. Bayley J-P, Oldenburg RA, Nuk J, Hoekstra AS, van der Meer CA, Korpershoek E, McGillivray B, Corssmit EPM, Dinjens WNM, de Krijger RR, Devilee P, Jansen JC, Hes FJ. Paraganglioma and pheochromocytoma upon maternal transmission of SDHD mutations. BMC Med Genet. 2014; 15:111. doi: 10.1186/s12881-014-0111-8.
49. Lefebvre M, Foulkes WD. Pheochromocytoma and paraganglioma syndromes: genetics and management update. Curr Oncol Tor Ont. 2014; 21:e8–17. doi: 10.3747/co.21.1579.
50. Feichtinger RG, Weis S, Mayr JA, Zimmermann F, Geilberger R, Sperl W, Kofler B. Alterations of oxidative phosphorylation complexes in astrocytomas. Glia. 2014; 62:514–25. doi: 10.1002/glia.22621.
51. Vidali S, Aminzadeh S, Lambert B, Rutherford T, Sperl W, Kofler B, Feichtinger RG. Mitochondria: The ketogenic diet—A metabolism-based therapy. Int J Biochem Cell Biol. 2015; 63:55–9. doi: 10.1016/j.biocel.2015.01.022.
52. Ju YS, Alexandrov LB, Gerstung M, Martincorena I, Nik-Zainal S, Ramakrishna M, Davies HR, Papaemmanuil E, Gundem G, Shlien A, Bolli N, Behjati S, Tarpey PS, et al. Origins and functional consequences of somatic mitochondrial DNA mutations in human cancer. eLife. 2014; 3. doi: 10.7554/eLife.02935.
53. Qing G, Skuli N, Mayes PA, Pawel B, Martinez D, Maris JM, Simon MC. Combinatorial regulation of neuroblastoma tumor progression by N-Myc and hypoxia inducible factor HIF-1alpha. Cancer Res. 2010; 70:10351–61. doi: 10.1158/0008-5472.CAN-10-0740.
54. Sausen M, Leary RJ, Jones S, Wu J, Reynolds CP, Liu X, Blackford A, Parmigiani G, Diaz LA, Papadopoulos N, Vogelstein B, Kinzler KW, Velculescu VE, et al. Integrated genomic analyses identify ARID1A and ARID1B alterations in the childhood cancer neuroblastoma. Nat Genet. 2013; 45:12–7. doi: 10.1038/ng.2493.
55. Pereira L, Soares P, Máximo V, Samuels DC. Somatic mitochondrial DNA mutations in cancer escape purifying selection and high pathogenicity mutations lead to the oncocytic phenotype: pathogenicity analysis of reported somatic mtDNA mutations in tumors. BMC Cancer. 2012; 12:53. doi: 10.1186/1471-2407-12-53.
56. Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, Clawson H, Spieth J, Hillier LW, Richards S, Weinstock GM, Wilson RK, Gibbs RA, et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes. Genome Res. 2005; 15:1034–50. doi: 10.1101/gr.3715005.
57. Pollard KS, Hubisz MJ, Rosenbloom KR, Siepel A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 2010; 20:110–21. doi: 10.1101/gr.097857.109.
58. Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003; 31:3406–15.
59. Picardi E, Pesole G. Mitochondrial genomes gleaned from human whole-exome sequencing. Nat Methods. 2012; 9:523–4. doi: 10.1038/nmeth.2029.
