The identification of mitochondrial DNA variants in glioblastoma multiforme

Mitochondrial DNA (mtDNA) resides exclusively within the matrix of the mitochondrion [1]. MtDNA is ~16.6 kb in size and encodes 13 proteins that are essential to the electron transfer chain, which produces the vast majority of cellular ATP through the biochemical process of oxidative phosphorylation (OXPHOS) [2]. It also consists of 22 tRNAs and 2 rRNAs and the D-Loop, which is the site of interaction for the nuclear-encoded factors that regulate mtDNA transcription and replication. The D-loop also contains two hypervariable regions, which are indicative of maternal ancestry [3].

The absence of intronic regions increases the likelihood of a variant residing within the coding regions of mtDNA [4], especially since the non-coding region accounts for only ~38% of the entire mitochondrial genome [5]. In the disease state, mutations are commonly observed throughout the mitochondrial genome [6, 7]. Due to the presence of multiple copies of mtDNA in the matrix of each mitochondrion and cells having large numbers of mitochondria, mutant and wild-type mtDNA can co-exist.

Primarily, aerobic glycolysis is central to tumor cell metabolism, as these cells exhibit three-times higher glycolytic activity than their normal counterparts [8]. This reflects the Warburg effect, whereby cancer cells undergo a metabolic transition, regardless of the oxygen concentration available [9]. This causes the cell to favour glycolysis rather than OXPHOS, as glycolysis generates ATP at a higher rate than OXPHOS, which produces ATP at a higher yield [10], and is thus more effective for cellular proliferation and growth. It is therefore likely that the presence of multiple mtDNA variants in a single cell or cohort of cells would promote tumorigenesis as there would be a decreased ability to generate ATP through OXPHOS and an increased reliance on glycolysis.

Generally, there is an accumulation in the number of mtDNA variants before a cell experiences detrimental or phenotypic effects [11] and this is likely to be the case for tumor cells. Nevertheless, it remains to be determined whether mtDNA variants present in tumors arise from the clonal expansion of existing variants or if they occur de novo. Normally, variants must be present at significant levels to allow for their detection [12]. However, next generation sequencing enables whole mitochondrial genomes to be sequenced at a high coverage of over 1000-fold, which provides sufficient depth to identify sequence variants at even low levels [13, 14].

Although the benefits of using such technology are clear, variants can be called that are interpreted as being true variants but are erroneous [15]. High resolution melt (HRM) analysis is an appropriate validation method to confirm next generation sequencing results by eliminating any errors that may have arisen from the sequencing reactions [16] or from long PCR amplification [17]. Confirmed variants can then be analysed by bioinformatics programs, such as SNPs & GO and MutPred that predict the impact of the SNPs on protein function [18].

Gliomas account for 80% of all primary malignant brain tumors and about 50% of gliomas give rise to Glioblastoma Multiforme (GBM) [19]. GBM has been rated most severe (level 4) on the World Health Organisation’s (WHO) scale of severity, due to its malignant nature [20]. For patients diagnosed with GBM, prognosis is poor with median survival estimated at 12–14 months [21]. Indeed, only about 3-5% of patients survive beyond 3 years [21].

A number of cell lines are representative of GBM. These include the frequently investigated lines, such as U87MG, SF-767 [22] and HSR-GBM1 [23]. Within patient GBM tissues, amplification of epidermal growth factor receptor (EGFR) genes is often observed [24]. This is evident in the GBM 6 and HK301 neurosphere cell lines, both of which express EGFRvIII. Others have been established from expanded cultures from high-grade human GBM tumors.

We have analysed a series of GBM cell lines using a combination of next generation sequencing and HRM analysis to identify variants within the mitochondrial genome that are indicative of tumorigenesis in GBM. We have then depleted one of the cell lines to varying degrees to determine whether the same variants reoccur following mtDNA replenishment as tumorigenesis progresses. We have applied several bioinformatic approaches to predict the severity of each mutation and to model the effects on protein structure. We have confirmed that these variants are indicative of GBM through their identification in GBM tumors. The frequency of variants is greatest within the D-loop and the origin of light strand replication in the non-coding regions whilst NADH dehydrogenase 6 (ND6) is the most susceptible gene for the development of mutations in the coding region.

This is the first study using next generation sequencing to identify the presence of mtDNA variants in GBM cell lines, which are representative of the relapse and/or refractory stage of the disease. We used a two-fold strategy to ensure that the variants identified were present within the cells and were not artifacts. Firstly, we amplified the two overlapping regions of the mitochondrial genome by long PCR followed by next generation sequencing, an approach used by others to, for example, identify the degree of heteroplasmy amongst different populations [29]. Importantly, this strategy ensured that we were not screening and, thus identifying variants, in nuclear pseudo-genes but rather targeting the mitochondrial genome itself [30]. Indeed, isolation of mitochondria and subsequent extraction of mtDNA can result in DNA samples possessing contaminating nuclear DNA. Secondly, we validated the presence of the variants using HRM on genomic DNA, which enabled us to confirm that the variants were ‘true’ calls, and not introduced during long PCR.

Overall, the GBM cell lines harbored proportionally more previously reported variants (16/29) than the normal, unaffected brain samples (20/60) following searches on the MitoMap variant database. Specifically, for the GBM cell lines, 4/13 variants that have not been previously reported resided within the non-coding region and 12/19 within the coding region. For the normal samples, there were two previously unidentified variants in the non-coding region and 18 in the coding region. Indeed, it was most apparent that there were only three variants common to the GBM cell lines and the normal samples, one in the D-loop at position 16519 bp (D-loop) and two in the coding genes at positions 14159 bp and 14160 bp (ND6). The 16519 bp variant has been identified in multiple studies [31] and was present in both cohorts at similar levels. However, the 14159 bp and 14160 bp variants have not been previously reported. Consequently, it appears that two distinct populations of variants define the two cohorts.

Although the variants identified in the GBM and normal cohorts appear to be randomly distributed throughout the mitochondrial genome, certain regions were not affected. No variants were identified in the NADH dehydrogenase 1 (ND1) gene or the 22 tRNAs suggesting that the incorporation of variants was highly selective to the other 12 coding genes and the D-loop region. Moreover, different regions of the mitochondrial genome were more susceptible to mutation in the normal samples than the GBM cell lines. For the GBM cells, ND6 and ND4 were the most susceptible genes in the coding region whilst the origin of light strand replication was the most susceptible non-coding gene. The ND6, ND4 and D-loop variants were also present in the majority of patient samples. For the normal samples, the D-loop was the most susceptible non-coding region whilst CYT B was the most susceptible gene in the coding region. Interestingly, for both cohorts, COX I was the least susceptible. COX I is one of the subunits of complex IV of the electron transfer chain protein and is amongst the most conserved of the proteins encoded by mtDNA [32]. Nevertheless, the COX II and Cytochrome C Oxidase III (COX III) genes, which also encode complex IV subunits, were more frequently susceptible to mutation than COX I in the GBM cell lines (2.26 fold) and the normal samples (7.87 fold), respectively. CYT B is also highly conserved [32] and the presence of variants in this gene would explain why they were predicted to be disease causing.

Although there was overlap between the various regions of the mitochondrial genome that were affected, there were regions that were specific to each cohort. The variants detected in the origin of light strand replication, COX II and ND4L were specific to the GBM cell lines, whereas variants specific to the normal samples were detected in the NADH dehydrogenase 2 (ND2), NADH dehydrogenase 3 (ND3), ATPase 6 (ATP6), ATPase 8 (ATP8) and COX III genes. This suggests that tumorigenesis targets specific regions of the mitochondrial genome for mutagenesis. To ascertain whether our outcomes were not biased by the cell lines and samples used, we confirmed through PHYLIP analysis that the cell lines and normal samples were from similar lineages and that their respective variants were not indicative of clustered cohorts. Furthermore, as many of the variants have not been previously reported, it appears that we have identified a series of mtDNA variants that are indicative of GBM, which can be validated in a large cohort.

Mitochondrial genetics is becoming increasingly important to our understanding of aging, metabolic disease, demographic patterns of human migration and cancer [33]. However, there is no association between mitochondrial dysfunction and the development of tumors in patients diagnosed with mitochondrial disease, where a single variant or deletion is present at high levels [34]. Moreover, no association has been effectively described for the existence of specific variants with a potential hotspot region for mutation [33]. Furthermore, others have argued that the nuclear genome is primarily responsible for the onset of tumorigenesis [35, 36]. However, multiple mtDNA variants scattered throughout the mitochondrial genome [37, 38] and in specific genes, such as ND5, are associated with the progression of tumorigenesis [39]. Nevertheless, studies using mouse transmitochondrial cybrids, where specific mtDNA variants have been introduced into cells with identical nuclear backgrounds, have shown that G→A at position 13997 bp and a frame-shift mutation at position 13885 bp caused by the insertion of a C residue, in the ND6 gene, were involved in metastasis and can affect tumor phenotype [40]. Other cybrid studies have also shown that mtDNA mutations act to positively induce tumorigenesis [41–45]. However, others hypothesise that there is active selection against unfavourable mtDNA mutations in tumor cells, due to the monoclonal origin of tumors whilst normal tissues are highly polyclonal in nature, thus promoting the selection of non-mutated mtDNA [46].

There are a number of studies that do not support the presence of mitochondrial mutations as a cause for tumorigenesis. In 2001, Kirches et al. analysed the D-loop of 17 paired GBM samples for mtDNA mutations, only to observe that 7 of the 17 tumors possessed variants [47]. The same authors later published a short report describing the absence of any sequence variant in 10 paired astrocytoma samples [48]. Moreover, in a study that focused on mtDNA mutations in mouse brain tumors, no pathogenic variants were detected across the mitochondrial genome in both induced and spontaneously developed tumors [49]. Seoane et al. also argued that the maintenance of mtDNA integrity is essential in order for tumor cells to survive increased reactive oxygen species production, which is associated with tumorigenesis [50, 51]. Altogether, it remains controversial whether mtDNA mutations arise as a cause or consequence of tumorigenesis. However, our data support the role of mtDNA as an initiator of tumorigenesis.

From our results, not only were the vast majority of HSR-GBM1 specific variants reestablished in cells depleted to varying degrees, but new variants were introduced as a function of mtDNA depletion. This is highlighted by those tumors generated from cells depleted to 20%, 3% and 0.2% of their original copy number. Interestingly, those variants that had been identified in the MitoMap database were present in all or nearly all of the newly formed tumors suggesting that there is a common core of mtDNA variants in tumors. Furthermore, 5/6 variants that had been established as de novo variants in the tumors have not been previously reported. To confirm that the newly derived variants were not indicative of clonal expansion from very few residual variants present in the cells, we analysed the raw data from the next generation sequencing runs for the HSR-GBM1 cell line and were able to confirm that these variants had not been called. Consequently, the induction of de novo mtDNA variants appears to be essential to the onset of tumorigenesis. Furthermore, our results reaffirm that ND4 and ND6 are the most susceptible genes to mutation in GBM and that GBM is likely to be a complex I disorder, which is supported by the one previous finding, prior to next generation sequencing, that complex I along with the D-loop is implicated in GBM [47]. Consequently, complex I would be a pharmaceutical target for the development of potential therapies.

Overall, the repopulation of depleted cells with mitochondria appears to be essential during tumorigenesis, as recently demonstrated [52]. The accumulation of mtDNA variants would result in a dysfunctional electron transfer chain that promotes reduced OXPHOS output, as is the case for GBM [53] and forces cells to adopt glycolysis, supported by low levels of OXPHOS. As a result, mtDNA copy number remains suppressed, which promotes the metabolic state of aerobic glycolysis, first proposed by Otto Warburg in 1956 [9]. This form of metabolism promotes self-renewal and cellular proliferation, which is also favoured by adult and embryonic stem cells [54–57]. However, embryonic and adult stem cells have the potential to differentiate into mature cell types, increase their mtDNA copy number and, consequently, switch from anaerobic metabolism to OXPHOS as their primary method of metabolism [52, 56, 58, 59].

We have identified a series of novel mtDNA variants by screening GBM cell lines. Mutations at positions 2130 bp within the 16 s rRNA and 5752 bp within the origin of light strand replication were sites of frequent mutation across the cell lines. Most of the GBM variants identified in the non-coding regions of mtDNA were within the D-loop region, whereas most variants in the coding portion resided within the ND4 region. However, ND6 was determined to be the most susceptible mtDNA region for the development of mutation in the GBM cells, given that it contained the most variants following normalization for the different sizes of the mtDNA regions. These outcomes were further endorsed by mtDNA depletion and in vivo replenishment during tumorigenesis and analysis of patient tumor samples, which support the hypothesis that mtDNA variants in complex I act as initiators of GBM.

KYY carried out the collection, assembly, analysis and interpretation of the data, and drafting of the manuscript; AD carried out the collection of data and analysis; JFD carried out the collection and analysis of all the in vivo data; GP designed and assisted in the analysis and interpretation of the protein modeling work; SJW carried out the collection of data; DKG carried out data analysis and provided financial support for the study; MM carried out the analysis and interpretation of the data; TGJ designed the in vivo work, carried out the analysis and interpretation of the data and provided the GBM cell lines and mice used in this study; JCSJ carried out the analysis and interpretation of the data, designed the experiments, wrote the manuscript and provided financial support. All authors read and approved the final manuscript.