Skip to main content


Modulation of mitochondrial DNA copy number in a model of glioblastoma induces changes to DNA methylation and gene expression of the nuclear genome in tumours



There are multiple copies of mitochondrial DNA (mtDNA) present in each cell type, and they are strictly regulated in a cell-specific manner by a group of nuclear-encoded mtDNA-specific replication factors. This strict regulation of mtDNA copy number is mediated by cell-specific DNA methylation of these replication factors. Glioblastoma multiforme, HSR-GBM1, cells are hyper-methylated and maintain low mtDNA copy number to support their tumorigenic status. We have previously shown that when HSR-GBM1 cells with 50% of their original mtDNA content were inoculated into mice, tumours grew more aggressively than non-depleted cells. However, when the cells possessed only 3% and 0.2% of their original mtDNA content, tumour formation was less frequent and the initiation of tumorigenesis was significantly delayed. Importantly, the process of tumorigenesis was dependent on mtDNA copy number being restored to pre-depletion levels.


By performing whole genome MeDIP-Seq and RNA-Seq on tumours generated from cells possessing 100%, 50%, 0.3% and 0.2% of their original mtDNA content, we determined that restoration of mtDNA copy number caused significant changes to both the nuclear methylome and its transcriptome for each tumour type. The affected genes were specifically associated with gene networks and pathways involving behaviour, nervous system development, cell differentiation and regulation of transcription and cellular processes. The mtDNA-specific replication factors were also modulated.


Our results highlight the bidirectional control of the nuclear and mitochondrial genomes through modulation of DNA methylation to control mtDNA copy number, which, in turn, modulates nuclear gene expression during tumorigenesis.


Mammalian cells have multiple copies of the mitochondrial genome (mtDNA) within each mitochondrion. The human mitochondrial genome is a circular, double-stranded genome that is 16.6 kb in size. It is essential for the production of cellular ATP, as it encodes 13 subunits of the electron transfer chain that conducts oxidative phosphorylation (OXPHOS) and is the cell’s major generator of ATP. It also encodes 22 transfer RNAs and 2 ribosomal RNAs [1].

mtDNA replication is driven by nuclear-encoded mtDNA-specific replication factors that translocate to the mitochondrion. The primary factor is the mtDNA-specific polymerase, polymerase gamma [2, 3], which is a heterotrimer enzyme composed of a catalytic subunit, subunit A (POLG), and two supporting subunits, subunits B (POLG2) [4]. The process of mtDNA replication is also supported by the mitochondrial helicase twinkle (TWNK), the single-stranded DNA binding protein (SSBP1) and the mitochondrial topoisomerase (TOP1MT) [5,6,7]. There are several other key factors including mitochondrial transcription factors A (TFAM), B1 (TFB1M) and B2 (TFB2M), which generate the precursor transcript used to initiate mtDNA replication [8, 9]. Furthermore, some upstream regulators of these mtDNA replication factors are highly involved in mitochondrial biogenesis and function. These include nuclear respiratory factors 1 and 2 (NRF1/2), peroxisome proliferator-activated receptor γ (PPARG) and its co-activator PGC1α (PPARGC1A), the Sirtuin family of genes (SIRT1-3), and the oestrogen-related receptors (ESRRA/B/G) [10,11,12,13,14].

mtDNA replication is strictly regulated during development and differentiation, which enables mature cells to acquire the requisite numbers of mtDNA copy to support their specific functions [15,16,17,18]. Initially, this is achieved by pluripotent (naïve) cells having established the mtDNA set point (reviewed in [19]). These naïve cells each possess around 200 copies of mtDNA, which promotes glycolysis as the favoured form of energy production, and cellular proliferation [20]. These copies of mtDNA serve as the initial template for mtDNA replication, allowing cells to acquire the appropriate numbers of mtDNA copy as they differentiate into mature cell types [16,17,18, 20, 21].

The inability to regulate mtDNA copy number affects cellular function and impedes developmental potential. For example, oocytes at the metaphase II stage that possess too few copies of mtDNA frequently fail to fertilize [22,23,24,25]. Likewise, certain types of cancer cells exhibit low mtDNA copy number and are unable to successfully complete differentiation [16, 18]. Indeed, cancer cells mainly rely on aerobic glycolysis for energy production, which allows for higher rates of cellular proliferation and prevents differentiation from taking place [16]. However, when they are induced to differentiate, they appear to be trapped in a ‘pseudo-differentiated’ state as they had not previously maintained the mtDNA set point and, as a result, the nuclear and mitochondrial genomes do not act in synchrony [19, 26]. Interestingly, when cancer cells are partially depleted of their mtDNA content, they undergo dedifferentiation and are then able to replicate mtDNA as they undergo differentiation, suggesting that they have re-established the mtDNA set point and synchrony between the nuclear and mitochondrial genomes [18].

During the early stages of development, the nuclear genome undergoes DNA demethylation mediated by the ten-eleven translocation methylcytosine dioxygenases (TET enzymes) to erase parental DNA methylation profiles. De novo DNA methylation then takes place to establish a new DNA methylation profile directed by the DNA methyltransferase family (DNMTs), which primes naïve cells for cellular differentiation [27, 28]. To this extent, cell-specific gene expression profiles are established during differentiation in synchrony with changes to DNA methylation patterns. Likewise, the DNA methylation status at exon 2 of POLG is a determinant of when this gene is expressed and, in turn, regulates mtDNA copy number in a cell-specific manner [15, 16]. This is supported by experiments using DNA demethylation agents, such as 5-azacytidine [29] and vitamin C [30], where modulation of DNA methylation at exon 2 of POLG increased mtDNA copy number in HSR-GBM1 cells derived from a glioblastoma multiforme (GBM) tumour [16, 18, 31].

The HSR-GBM1 cell line is a high-grade malignant GBM cell line that is characterized as being similar to stem-like neural precursors and is extensively DNA methylated, which contributes to its tumorigenic gene profile [32, 33]. However, its hyper-methylated profile is not established by the overexpression of the isocitrate dehydrogenases (IDH1/2) that harbour onco-mutations, as the alleles for these genes are wild type [34, 35]. Under normal circumstances, IDH enzymes act on the citric acid cycle to generate α-ketoglutarate, which is a co-factor of the TET enzymes that modulate DNA demethylation patterns [36,37,38,39,40]. However, overexpression of and mutations to the IDH genes in GBM result in a metabolic switch that produces 2-hydroxyglutarate and restricts DNA demethylation induced by the TET enzymes [36,37,38,39,40]. Consequently, HSR-GBM1 cells enable the analysis of modifications to DNA methylation profiles to be undertaken whereby the DNA methylation status of the cells is not influenced by mutations to key regulators of DNA demethylation and thus allows the effects of mtDNA copy number to be studied independently of these influences.

Interestingly, mtDNA depletion of HSR-GBM1 cells to varying amounts of mtDNA copy number affected tumour progression and frequency when these cells were inoculated into mice [18]. Progression and frequency were greatest in cells depleted to 50% of their original content, but tumour formation was less frequent and took significantly longer when cells possessed only 3% and 0.2% of their original mtDNA content [18]. Notably, mtDNA copy was restored to similar levels during in vivo tumorigenesis accompanied by DNA demethylation at exon 2 of POLG [17].

In order to determine whether global DNA methylation profiles were modulated following the restoration and maintenance of mtDNA copy number in end point tumours, we investigated the DNA methylation profiles of GBM tumours derived from HSR-GBM1 cells that possessed varying amounts of mtDNA copy number and exhibited different frequencies and progression in tumour formation. We used whole genome methylated DNA immunoprecipitation (MeDIP)-Seq. We matched the modulated regions with their transcriptional profiles to focus on their effects on gene expression. We also investigated the mtDNA replication factors to determine how they responded to the newly established interactions between the nucleus and the mitochondrial genome. Our results highlight the bidirectional control of the nuclear and mitochondrial genomes through modulation of DNA methylation to control mtDNA copy number and gene expression in tumours using the HSR-GBM1 cell line as a model.


Replenishment of mtDNA copy number

Tumours were previously generated from HSR-GBM1 cells possessing different levels of mtDNA copy number, namely GBM100 (possessing 100% of their mtDNA content), GBM50 (50% mtDNA content), GBM3 (3% mtDNA content) and GBM0.2 (0.2% mtDNA content) cells [18]. Tumours arising from GBM3 and GBM0.2 cells had significantly delayed initiation and reached end point at 83 and 90 days, respectively. They had a lower frequency of tumour formation (6/12 and 2/12, respectively) compared with GBM100 tumours, which took 65 days to reach end point with tumours forming from 11/12 inoculations [18]. GBM50 tumours had accelerated formation (61 days) [18]. Furthermore, mtDNA copy number was replenished in each of the end point tumours to pre-depletion levels without significant differences (Fig. 1), which highlights the need for sufficient mtDNA copy number to promote tumorigenesis, and likely explains why delayed tumour formation was observed in cells with lower levels of mtDNA copy number.

Fig. 1

Replenishment of mtDNA copy number in the GBM50, GBM3 and GBM0.2 tumours. Fold change in mtDNA copy number for GBM50, GBM3 and GBM0.2 tumours relative to the GBM100 tumours. GBM100, GBM50, GBM3 and GBM0.2 tumours were harvested at 65 days, 61 days, 83 days and 90 days after inoculation, respectively. Statistical significance was determined using one-way ANOVA (n = 3 tumours per group; technical triplicates per sample). No significant differences were identified between the groups

Modifications to DNA methylation resulting from re-establishment of mtDNA copy number during tumorigenesis

In order to determine whether DNA methylation of the nuclear genome was modulated following re-establishment of mtDNA copy number during tumour progression, we assessed the global DNA methylation profiles of each of the tumour types at end point. Whole genome MeDIP-Seq was performed on DNA extracted from individual tumours from each cohort to identify key regions that might be modulated. Using the MEDIPS package, the quality of the sequencing reads was determined using the saturation and enrichment analyses and each achieved high standards, indicating that there was sufficient coverage and enrichment of the methylated DNA (Additional file 1). Regions with total coverage of 100 × (minrowsum = 100) or greater reads were used for statistical analysis, which is deemed to be a more stringent filter than the default threshold of tenfold coverage [31].

A common feature of cancer cells is genomic instability, which causes significant copy number variation (CNV) when compared to healthy controls [41]. CNV can also affect the outcomes from comparative analysis after MeDIP-Seq, which, if ignored, can introduce false positive results due to CNVs inducing higher or lower levels of enrichment rather than being reflective of the real differences in DNA methylation [41]. Using the Illumina Global Screening Array, no CNVs were identified to be significantly different between the tumour groups, as assessed using the NEXUS Copy Number Module. This negates the potential effect of CNV on the analysis of MeDIP-Seq outcomes (see Additional file 2).

Differentially methylated regions (DMRs) were identified between the GBM100 tumours and tumours formed from cells possessing varying levels of mtDNA. Table 1 shows the number of DMRs identified in the GBM50, GBM3 and GBM0.2 tumours when compared with the GBM100 tumours using different filters for statistical significance. With the increase in statistical power (p value) from ≤ 0.05, 0.01 to 0.001, fewer DMRs were identified in each comparison. Overall, there were fewer DMRs identified in the GBM3 tumours followed by the GBM50 tumours whilst there were far more DMRs identified in the GBM0.2 tumours. Using a cut-off p value of 0.001, for instance, 124 DMRs were identified between the GBM50 tumours and the GBM100 tumours, 36 DMRs for the GBM3 tumours and 16,565 DMRs for the GBM0.2 tumours. The greater number of DMRs identified for the GBM0.2 tumours indicated that the DNA methylation profile for this cohort of tumours was extensively modified as a result of the increased levels of mtDNA depletion in their initiating cells and subsequent restoration during tumorigenesis. We further applied the Bonferroni correction to adjust p values to control for false positive calls. After correction, there were 143 DMRs identified in the GBM0.2 tumours, but there were no significant DMRs identified in the GBM50 and GBM3 cohorts of tumours.

Table 1 The number of DMRs identified in tumours (n = 3 per group) formed from depleted cells compared with non-depleted cells

Both hyper- and hypo-methylated states were observed for tumours derived from cells with lower levels of mtDNA copy number. However, the majority of the DMRs were hypo-methylated compared to the GBM100 tumours (Table 1). There were 31 hyper-methylated DMRs (21.68%) identified in the comparison between GBM0.2 and GBM100 tumours using an adjusted p value of 0.05. Indeed, the tendency to hypo-methylation was also observed in the other tumours. Taking the results identified using a p value of, for instance, 0.001, only 7.26% of the DMRs in the GBM50 tumours and 44.44% of the DMRs in the GBM3 tumours were hyper-methylated. The different proportions of hypo- and hyper-methylation further suggest that when tumorigenesis was initiated from different levels of mtDNA copy number, different patterns of DNA methylation were induced.

The annotation of the DMRs and gene ontology analysis for the DMR-overlapping genes

Firstly, the annotation of genomic regions was performed on the four groups of DMRs identified using the filters of a p value ≤ 0.001 and an adjusted p value ≤ 0.05 [GBM50 vs GBM100 (p ≤ 0.001), GBM3 vs GBM100 (p ≤ 0.001), GBM0.2 vs GBM100 (p ≤ 0.001) and GBM0.2 vs GBM100 (adjusted p ≤ 0.05)]. Overall, the results showed that over 49% of the DMRs were enriched in intragenic regions, namely exons, introns and promoter regions, and they were the most enriched regions (Fig. 2a–d). Each of these four groups of DMRs showed the greatest enrichment at introns, which was 52.51% of the GBM50 DMRs (Fig. 2a), 39.22% of the GBM3 DMRs (Fig. 2b), 34.73% of the GBM0.2 DMRs (Fig. 2c; p ≤ 0.001) and 35.94% of the GBM0.2 DMRs (Fig. 2d; adjusted p ≤ 0.05). For the GBM50 and GBM3 DMRs, the second most enriched regions were the distal intergenic regions, which were enriched at 22.35% and 37.25%, respectively (Fig. 2a, b). Interestingly, for the GBM0.2 DMRs, identified using the filters of p ≤ 0.001 and adjusted p ≤ 0.05, similar patterns were observed regardless of the numbers of DMRs (Fig. 2c, d). In contrast to the GBM50 and GBM3 DMRs, the second most enriched region for the GBM0.2 DMRs was the promoter regions, which were 23.60% and 25.81% of the DMRs after filtering at p ≤ 0.001 and adjusted p ≤ 0.05, respectively (Fig. 2c, d). The distal intergenic regions were ranked third accounting for 16.51% and 18.89% of the GBM0.2 DMRs identified using the filters of p ≤ 0.001 and adjusted p ≤ 0.05, respectively (Fig. 2c, d). Untranslated regions (UTRs), downstream regions of genes and super enhancers accounted for less than 14% of the DMRs in each comparison. Interestingly, over 80% of the identified DMRs overlapped with topologically associating domains (TADs) that interact with other genomic regions and are considered to be the regulatory units of the genome [42].

Fig. 2

Annotations of the DMRs. Based on their genomic locations, DMRs identified in the comparison: a GBM50 versus GBM100 (p ≤ 0.001); b GBM3 versus GBM100 (p ≤ 0.001); c GBM0.2 versus GBM100 (p ≤ 0.001); and d GBM0.2 versus GBM100 (adjusted. p ≤ 0.05) were classified into promoters, exons, introns, UTRs, downstream regions, distal intergenic regions and super enhancers using the ChIPSeeker package (n = 3 tumours per group). DMR-overlapping regions in brown signify promoter regions, violet-red indicates exons, green indicates introns, blue indicates promoters, teal represents UTRs, blue represents downstream regions, orange indicates intergenic regions, and pink indicates super enhancers. The percentage of each group is indicated in the key to each figure

To understand the related functions of the DMRs, functional classification and statistical overrepresentation analyses were performed on the DMR-overlapping genes using the PANTHER Classification System. The functional classification of the DMR-overlapping genes showed that the GBM50 and GBM3 tumours mostly clustered into cellular and metabolic processes (Fig. 3). Furthermore, the pathways of biological regulation, localization, response to stimuli and developmental processes were present in both cohorts (Fig. 3). Due to the large number of DMR-overlapping genes identified in the GBM0.2 tumours using the two different filters, statistical overrepresentation analysis was, therefore, performed. Interestingly, the biological processes enriched in the DMR-overlapping genes showed strong associations with developmental processes, especially in the DMR-overlapping genes identified with an adjusted p value of 0.05 (Table 2). The common developmental pathways identified in both analyses (shown in italics) included developmental, nervous system development and cell differentiation processes (Table 2). This is in line with a previous finding that, after partial mtDNA depletion with 7 and 14-day ddC treatment, differentiation was able to take place in the GBM cells [18]. Moreover, regulation of transcription from the RNA polymerase II promoter was found to be commonly modulated, which reflected the potential for transcriptional changes that could result from the DMRs.

Fig. 3

Functional classification analysis of the DMR-overlapping genes in the GBM50 and GBM3 tumours. DMR-overlapping genes in: a GBM50 tumours; and b GBM3 tumours were classified into subgroups based on their associated biological processes using the PANTHER system

Table 2 Pathways affected by the DMR-overlapping genes in the GBM0.2 tumours using the PANTHER Classification system

Validation of the gene expression of the regulators of DNA methylation

DNA methyltransferases (DNMTs), namely DNMT1/3A/3B, are responsible for catalysing DNA methylation at cytosines by converting them to 5-methylcytosines (5mC). On the other hand, the ten-eleven translocation (TET) enzymes, namely TET1/2/3, can demethylate 5mC to 5-hydroxymethylcytosine (5hmC) [27, 28]. The changes to these DNA methylation and demethylation factors can provide insights into how the changes to the DNA methylation profiles occurred. Therefore, the gene expression levels of the DNA methylation regulators were assessed using a custom-made Fluidigm qPCR array (Fig. 4). There was ≥ 1.5-fold change in expression of DNMT1 in GBM50 and GBM0.2 tumours, suggesting that the reduction in mtDNA copy number strongly induced DNA demethylation in the tumours, and is probably the major reason why the majority of regions were hypo-methylated (Table 1). DNMT3A was also significantly down-regulated by over 1.5-fold in the GBM50 tumours (Fig. 4). None of the DNA methylation factors was found to be differentially expressed in the GBM3 tumours compared with the GBM100 tumours, which further supported the results that fewer DMRs were identified in this group. The expression of TET1 showed a trend of up-regulation only in the GBM0.2 tumours. As there was a higher percentage of hyper-DMRs in the GBM0.2 tumours, the up-regulation of TET1 might be a response to counter the increased hyper-methylation observed.

Fig. 4

Differential expression of the regulators of DNA methylation. Bars represent the mean of the relative quantification levels relative to the GBM100 tumours (relative expression = 1). Error bars show SEM (n = 3 tumours per group; technical replicates per sample). Statistical significance was determined by One-way ANOVA. *, ** and *** p values < 0.05, 0.01 and 0.001, respectively

Furthermore, it has been reported that the genomes of low grade and secondary glioblastomas are extensively methylated to support the tumorigenic transcriptome that they exhibit [32, 38,39,40, 43]. This is further supported by the findings that the overexpression of the mutated isocitrate dehydrogenase (IDH1/2) in the citric acid cycle inhibits the activity of TET2 and thus contributes to the hyper-methylated genome [36,37,38,39,40]. However, the HSR-GBM cell line possesses wild-type alleles for IDH1/2 [34, 35], which makes it an excellent model to investigate the impact of mtDNA content in tumour-initiating cells on DNA methylation patterns during tumorigenesis, especially as non-mutated IDH1 has been shown to promote cell growth in GBM [44]. We, therefore, investigated the expression of IDH1/2 in these tumours and found that expression for both genes was down-regulated in the GBM50 tumours (Fig. 4), which likely reduced their ability to induce hyper-methylation and supports the hypo-methylation observed in these tumours.

DMR-overlapping differentially expressed genes and differentiation markers after restoration of tumorigenic capacity

To investigate which genes changed their patterns of expression in response to the altered levels of DNA methylation, we compared the identified DMRs with the differentially expressed genes previously reported on the same samples [17]. We focused on the DMRs identified using an adjusted p value of 0.05, which were considered to be the genomic regions that had undergone most significant changes in DNA methylation. There were 9 DMR-overlapping differentially expressed genes, namely BAIAP2, L3MBT1L, KCNC1, GPSM1, SLC27A1, MAF, OGFR, MICALL2 and RHOT2 (Fig. 5a). The changes to the levels of gene expression were further validated for these genes (Fig. 5a). Interestingly, the GBM3 tumours, with the least number of DMRs, showed more variable patterns in expression amongst the DMR-overlapping genes when compared with the GBM100 tumours. RHOT2 and OGFR were significantly up-regulated in the GBM3 tumours (p < 0.05 and p < 0.01, respectively). The greater modification that occurred in the DNA methylation profiles of the GBM50 and GBM0.2 tumours promoted gene expression profiles more similar to the tumorigenic features of the GBM100 tumours, whereas less modification in the GBM3 tumours failed to mimic the transcriptional changes in these DMR-overlapping genes.

Fig. 5

Transcriptional profiles of DMR-overlapping differentially expressed genes and markers of neurogenesis and neural stem cells. a Differential expression of the DMR-overlapping differentially expressed genes. Bars represent the mean of the relative quantification levels relative to the GBM100 tumours (relative expression = 1). Statistical significance was determined by One-way ANOVA. Error bars show SEM (n = 3 per group; technical triplicates per sample). *, ** p values < 0.05 and 0.01, respectively. b Heatmap of markers of neurogenesis and neural stem cells from the RNA-Seq profiles. For the GBM50, GBM3 and GBM0.2 cohorts, the fold changes (log2) of read counts to the mean value of the GBM100 cohort were plotted. The colour scheme from blue, white to red represents the level of expression from twofold down-regulation to twofold up-regulation

As the DMRs were found to be tightly associated with developmental processes (Fig. 3 and Table 2), we also focused on the neural markers and genes in neurogenesis that have been reported in cells having undergone mtDNA depletion to the same levels as those used to form tumours [18]. Indeed, it has been previously shown that the partial mtDNA depletion of GBM cells could lead to increased expression of genes associated with early developmental processes [18]. In tumours formed from these cells, an overview of gene expression identified by RNA-Seq data is shown (Fig. 5b). Both up-regulation and down-regulation were observed. To make a comparison between the in vivo tumours and in vitro cells, the negative regulator of cell proliferation ALK was found to be down-regulated in the mtDNA-depleted tumours. However, ALK was found to be up-regulated in the mtDNA-depleted cultured cells. Likewise, the regulators of cell differentiation and proliferation, VEGFA, DLG4, CDK, NRP2 and ACHE, were also found to be down-regulated in the mtDNA-depleted tumours but presented higher levels of expression in the cells. The regulators of transcription and differentiation, ASCL1 and DLL1 in the Notch signalling pathway, were found to be up-regulated in the mtDNA-depleted tumours but down-regulated in the cells. The regulator of synaptogenesis, APOE, was found to be up-regulated in the GBM50 and GBM3 tumours, but down-regulated in the GBM0.2 tumours and the mtDNA-depleted cells. Moreover, FGF13, SEMA4D, PAPD6B and NRCAM were found to be commonly up-regulated in both cells and tumours. As the cells recovered their tumorigenicity after restoration of mtDNA copy number in vivo, the shifts in gene expression between the tumours and the cells indicated that the restoration and maintenance of mtDNA copy number in tumorigenesis is tightly associated with the regulation of differentiation.

Modulation of the mtDNA replication factors for the maintenance of mtDNA copy number in tumorigenesis

As mtDNA copy number was replenished to similar levels in each of the tumour types (Fig. 1), we determined how the mtDNA-specific replication factors were affected. In total, 20 known mtDNA replication factors were screened for differentially methylated intragenic CpG islands (CGIs). The fold changes in relative methylation scores for the GBM100 tumours of 28 promoter and gene-body CGIs within the mtDNA replication factors were plotted (Fig. 6a). An increasing trend in fold changes was observed from the GBM50 to the GBM0.2 tumours. Particularly, for the CGIs, higher levels of methylation were found in the GBM50 tumours compared to the GBM100 tumours (shown in red blocks) and the fold changes generally increased to an even higher level in the GBM0.2 tumours. This trend was mostly observed amongst the gene-body CGIs. For the CGIs that were more hypo-methylated in the GBM50 tumours than the GBM100 tumours (shown in blue blocks), the levels of DNA demethylation were indicative of minor shifts in the GBM0.2 tumours.

Fig. 6

Overview of the changes to DNA methylation and transcription of the mtDNA replication factors. a DNA methylation levels of the intragenic CGIs associated with the mtDNA replication factors. Fold changes in relative methylation scores to the GBM100 tumours were plotted for the three biological replicates from the GBM50, GBM3 and GBM0.2 cohorts. The colour scheme from blue, white to red represents the level of DNA methylation from low to high. “P”, “E” and “I” following each gene name indicates CGIs that are located at the promoter regions, exons and introns, respectively. b Differential expression of the mtDNA replication factors. Bars represent the mean of the relative quantification levels relative to the GBM100 tumours (relative expression = 1). Error bars show SEM (n = 3 tumours per group; technical replicates per sample). Statistical significance was determined by one-way ANOVA. *, ** and *** p values < 0.05, 0.01 and 0.001, respectively

Furthermore, the gene expression levels were also assessed using a Fluidigm qPCR array (Fig. 6b). Overall, the levels of expression of the mtDNA replication factors in the tumours formed from cells possessing varying levels of mtDNA were found to be down-regulated, which likely restricts mtDNA replication during tumorigenesis after the restoration of mtDNA copy number to its original levels. ESRRA, NRF1, POLG, PPARGCA1 and SIRT1 were identified to be significantly down-regulated in all of the mtDNA-depleted tumours. The mtDNA replication factors POLG2, TOP1 MT and TWNK were also found to be down-regulated in the GBM50 and GBM3 tumours. Moreover, the key mtDNA transcription factors POLRMT, TFAM, TFB1M and TFB2M were down-regulated in the GBM50 and GBM3 tumours. NRF2 and SIRT3 were only down-regulated in the GBM50 tumours. The hypoxia regulator STAT3 was found to be significantly down-regulated in the GBM50 and GBM0.2 tumours. ESRRG was only down-regulated in the GBM3 tumours. This indicates that mtDNA replenishment had been completed and the low levels of expression of the mtDNA replication factors were indicative of mtDNA turnover during cell division rather than the active repopulation of depleted cells with mtDNA.


We have shown that bidirectional cooperation between the nuclear and mitochondrial genomes is important for tumour formation. To this extent, the DNA methylation and gene expression profiles of the nuclear genome were extensively modulated following the restoration of mtDNA copy number to pre-depletion levels in tumours that developed to end point. Whilst it is conceivable that the mtDNA depletion process did not reduce mtDNA copy in individual cells from each group uniformly, it is unlikely that this occurred. If a small population of cells with higher mtDNA copy number had been selected for, there would have been greater similarity amongst the DNA methylation and gene expression profiles across the groups. However, cells from each group exhibited distinct gene expression profiles prior to inoculation [18] and DNA methylation and gene expression profiles as end point tumours.

The restoration of mtDNA copy number can be achieved through the cell’s own mtDNA replication machinery or the horizontal transfer of mitochondria from surrounding cells [17, 45, 46]. It was previously shown that the same tumours possessed only human mtDNA and not mouse mtDNA from the surrounding stroma [47]. However, the restoration of mtDNA copy number also involved the accumulation of mtDNA variants, which increased in number as a function of the degree of restoration of mtDNA copy number. Each of the tumours exhibited a gain in de novo variants primarily associated with ND4 and ND6, which encode subunits of complex I of the electron transfer chain, and the D-loop. However, these variants have been identified in other GBM cell lines [47]. Consequently, the presence and modulation of mtDNA are important to tumorigenesis given that cells require mtDNA to be fully functional [48].

The effects of DNA methylation on transcription vary dependent on the location of the methylated sites. Methylated CGIs within the promoter regions are known to repress transcription [49], whereas methylation within gene bodies correlates positively and negatively with transcriptional elongation [50]. In our work, the degree of restoration of mtDNA copy number differentially affected the levels of DNA methylation in the end point tumours. However, the majority of the DMRs were hypo-methylated in each comparison. Whilst there was high enrichment in the intragenic regions, namely promoters, exons and introns, there were differences in the numbers of DMRs affected and their distribution across the nuclear genome, especially for the promoter and exon regions amongst the comparisons. The DMRs in the GBM50 and GBM0.2 tumours also overlapped with several enhancers, which are documented to contribute to tumour progression and cancer cell plasticity [51]. Furthermore, over 80% of the DMRs overlapped with TADs, the regulatory units of the genome [42], which suggests that their modulation could affect the structural formation and, therefore, genetic activities of the genome [42].

These outcomes are likely mediated by the down-regulation in levels of expression of the DNA methyltransferases, especially DNMT1. For the GBM0.2 tumours, DNA demethylation was further enhanced by the up-regulation of TET1, which mediates the transition from 5mC to 5hmC [52], resulting in these tumours being most hypo-methylated. IDH1 and IDH2, the cytoplasmic and mitochondrial isoforms respectively, had lower levels of expression in the tumours derived from mtDNA-depleted cells than GBM100 tumours. This suggests that α-ketoglutarate could continue to act as a co-factor of the TET enzymes and did not induce a switch to 2-hydroxyglutarate to promote hyper-methylation [43, 44, 53]. This is consistent with HSR-GBM1 cells possessing wild type copies of the IDH genes [34, 35], which is common in high-grade (IV) astrocytomas [38, 40, 44]. One-carbon metabolism also takes place in the mitochondria resulting in the generation of S-adenosylmethionine, the universal methyl group donor [54]. S-adenosylmethionine is present at abnormal levels in cancer cells and could affect DNA methylation profiles during tumorigenesis [54]. We would anticipate that S-adenosylmethionine activity was down-regulated as a result of mtDNA depletion, which reduces mitochondrial function [55], and, thus, promotes the hypo-methylated state of the resultant tumours.

During early development, large-scale DNA demethylation takes place to reset the nuclear genome to a naïve state [27, 28]. This is coupled with the strict regulation of mtDNA copy number to establish the mtDNA set point [15]. As a result, the two genomes act in tandem to overcome a number of molecular check points at different stages of development, which include mtDNA replication turnover events [56]. Indeed, Polg-/- mice do not possess sufficient mtDNA at E6.5 and cannot initiate mtDNA replication at this key developmental check point and, consequently, die [57]. Likewise, somatic cells reprogrammed to a naïve, pluripotent state, similar to embryonic stem cells, can fail to re-establish the mtDNA set point and lose their differentiation potential [56]. HSR-GBM1 cells are cancer stem cells that express key early neural markers, maintain low mtDNA copy number and use glycolysis for energy production, which are typical facets of a stem cell [33]. However, when they are induced to differentiate into high OXPHOS-derived ATP requiring cells, for example neurons and astrocytes, they fail to replicate mtDNA in synchrony with changes to the nuclear genome and their progress stalls [56]. Nevertheless, mtDNA partially depleted HSR-GBM1 cells can re-establish their mtDNA set point and expand their mtDNA copy number in order to complete differentiation [18] (please refer to Fig. 7 ‘in vitro’). However, the degree of restoration of mtDNA copy number influenced the differentiation and transcriptional potentials of the cells as determined by RNA-Seq analysis of end point tumours. The most affected networks between the GBM100 and the GBM0.2 and GBM3 tumours were cellular development, cellular growth and proliferation, whilst cancer, cell cycle and cellular development were most affected in the GBM50 tumours [17].

Fig. 7

Modulation of mtDNA copy number affects the fate of HSR-GBM1 cells. In the in vitro model of HSR-GBM1 cells, depletion of mtDNA (yellow arrow) for 7/14 days results in demethylation of the nuclear genome (green) and rescues differentiation to astrocytes that is otherwise blocked (green arrow). In the in vivo model of HSR-GBM1 cells, varying levels of depletion result in different frequencies of progression to tumorigenesis (orange arrow). 50% depletion accelerated tumorigenesis (thicker orange arrow), whereas further depletion slowed the initiation of tumorigenesis by modulating differentiation networks, nuclear gene expression and the mtDNA replication factors

Although the GBM50 tumours are more aggressive, they appear to restrict mtDNA replenishment to pre-depletion levels, which ensures that the balance between mtDNA copy number and the methylation of the chromosomal genome is at levels that promote tumorigenesis, as is the case for the GBM3 and GBM0.2 tumours. From a mtDNA perspective, there is increased hyper-methylation in the majority of gene-body CGIs of the mtDNA replication factors of the mtDNA-depleted tumours compared with the GBM100 tumours. This is exemplified by increased methylation at exon 2 of POLG and exon 8 of TOP1MT and decreased levels of expression for these two genes, which are known to regulate mtDNA copy number in a cell-specific manner through DNA methylation at these intragenic regions [15, 16, 31]. Interestingly, the use of DNA methylation agents, such as 5-Azacytidine and Vitamin C, on cultured HSR-GBM1 and reprogrammed somatic cells can restore mtDNA replication turnover events and promote differentiation [16], which likely explains the use of these agents in clinical settings for cancer patients [58, 59].

The DMR-overlapping genes affected a variety of biological processes. In the GBM50 and GBM3 tumours, the primary effects were on cellular and metabolic processes. In the GBM0.2 tumours, developmental processes were more overrepresented, including cell differentiation and nervous system development. Of the nine DMR-overlapping differentially expressed genes, four were specific to neurogenesis. SLC27A1 was significantly down-regulated in a SOX2-knockdown GBM cell line characterized by abolished dedifferentiation and decreased tumorigenesis [60]. Mutations in KCNC1 have been reported as the driving cause for progressive myoclonus epilepsies [61]. Indeed, a number of DNA methylation signatures have been identified in several neurodevelopmental syndromes, including Coffin–Siris, Kabuki and CHARGE syndromes, which are implicated in a variety of cancers [62, 63].

Although the GBM3 tumours presented with the least number of DMRs, they had more variable patterns of expression amongst the DMR-overlapping genes. For example, RHOT2 and OGFR were significantly up-regulated in the GBM3 tumours. RHOT2 is a miro GTPase that regulates mitochondrial transport, distribution and dynamics, especially in neurodegenerative disorders [64, 65], and is differentially expressed in different types (CD133+ and CD133) of GBM stem cells [66]. It is usually only expressed in type I GBM cancer stem cells and is co-expressed with the stem cell markers SOX2, SOX11 and OLIG2 [66]. This indicates the expression of RHOT2 is important for regulating mitochondria and their dynamics during early development, which requires the presence of mtDNA given that cells need mtDNA to be fully functional [48]. On the other hand, OGFR is a negative regulator of cell proliferation and reduces proliferation of astrocytes in cell culture [67]. Consequently, GBM3 tumours differentially expressed two genes associated with different stages of development unlike the other tumour types. Nevertheless, the greater levels of modification that occurred in the DNA methylation profiles of the GBM50 and GBM0.2 tumours promoted gene expression profiles more indicative of the tumorigenic features of the GBM100 tumours. A similar situation has been observed in tumour cell lines possessing the same chromosomal background but different mtDNA haplotypes where early tumours from each haplotype expressed a distinct set of genes [68]. However, in this case, a series of commonly expressed genes were expressed again highlighting how the mitochondrial genome can affect tumorigenesis.

It is evident that GBM0.2 cells, which had undergone a longer period of mtDNA depletion through an agent that specifically targets POLG from interacting with the mitochondrial genome [69], adopted a more dedifferentiated state than the partially depleted cell types, as demonstrated by the expression of early neural developmental genes [18]. Indeed, many key regulators of cell proliferation and differentiation were down-regulated in the mtDNA-depleted tumours, but up-regulated in the mtDNA-depleted cells [18], whilst regulators of transcription and differentiation in the Notch signalling pathway were up-regulated in the tumours derived from mtDNA-depleted cells but down-regulated in the cells [18]. Likewise, the regulator of synaptogenesis, APOE, was up-regulated in the GBM50 and GBM3 tumours but down-regulated in the mtDNA-depleted cells [18]. Consequently, since the mtDNA-depleted cells recovered their tumorigenicity as mtDNA copy number was restored in vivo, it is evident that the restoration and maintenance of mtDNA copy number in tumour-initiating cells is tightly associated with the regulation of neural differentiation in glioblastoma.

The regulation of mtDNA copy number, DNA methylation and nuclear gene expression at different stages of development highlights the need for cells to establish a balance between these three components to promote their specific fates. In the case of tumours derived from cells with different levels of mtDNA, unique DNA methylation profiles are established that regulate gene expression and mtDNA copy number in end point tumours (see Fig. 7 ‘in vivo’). To this extent, GBM0.2 cells attempted to modulate global DNA methylation patterns and transcription of the mtDNA replication factors and differentiation markers at the expense of tumorigenesis resulting in fewer tumours forming and requiring significantly longer to form. However, GBM50 cells were more tumorigenic, as they had struck the appropriate balance between mtDNA copy number and the methylated state of the chromosomal genome to promote tumorigenesis. The proposed model shown in Fig. 7 indicates how modulation of mtDNA copy number in tumour-initiating cells induces changes to DNA methylation of the nuclear genome and, therefore, can affect the fate of tumour cells in in vitro and in vivo environments. This highlights the synergy required between the two genomes in establishing tumorigenesis.


In all, we have shown that the restoration of mtDNA copy number during tumorigenesis induces major changes to the nuclear genome that resulted in differential DNA methylation and expression of genes. These changes enriched developmental processes and essential metabolic pathways associated with GBM. In addition, the changes to the nuclear-encoded mtDNA replication factors highlight the synergy between the nuclear and mitochondrial genomes in restoring tumorigenic capacity. This was clearly demonstrated in the GBM0.2 tumours formed from cells having undergone the highest levels of depletion and requiring the longest to initiate tumour formation. They underwent more extensive DNA methylation remodelling at key CGIs and within the intragenic regions of the mtDNA replication factors to maintain similar transcriptional levels. Whilst we recognize that we have only focused on the HSR-GBM1 tumour model, our results highlight the bidirectional control of the nuclear and mitochondrial genomes through modulation of DNA methylation in response to mtDNA copy number and to control gene expression in tumorigenesis.


Cell culture and Xenograft models

GBM tumours were previously generated from HSR-GBM1 cells possessing different levels of mtDNA copy number after treatment with 10 μm 2′–3′-dideoxycytidine (ddC), a mtDNA depletion agent that directly inhibits the interaction of POLG with the mitochondrial genome [69], in the presence of 50 mg/mL uridine (Sigma-Aldrich, MO, USA). Cells were depleted to 50% (GBM50), 3% (GBM3), 0.2% (GBM0.2), and 100% (GBM100) of their original mtDNA content, as described in [18]. To this extent, HSR-GBM1 cells were cultured in complete neural stem cell media consisting of Dulbecco’s Modified Eagle Medium/Nutrient Mixture (DMEM/F-12) Media (Thermo Fisher Scientific, MA, USA), 2% StemPro neural supplement (Thermo Fisher Scientific), 20 ng/mL basic fibroblast growth factor (bFGF; Merck Millipore MO, USA) and 20 ng/mL epidermal growth factor (EGF; Merck Millipore) at 37 °C, 5% CO2 and 95% humidity.

The animal work was approved by the Animal Ethics Committee, Monash University, Approval Number: MMCA/2011/76. Briefly, 0.5 million HSR-GBM1 tumour cells in 100 mL of medium were inoculated subcutaneously into both flanks of 5- to 6-week-old, female BALB/c nude mice (Animal Research Centre, Perth, Australia). Tumour growth rates and volumes were reported in [18]. For each group, cells were injected into 12 mice to form tumours, as detailed in [18]. In all, 11 tumours formed from GBM100 cells, 10 from GBM50 cells, 6 from GBM3 cells and 2 from GBM0.2 cells [18]. GBM100 tumours reached an average volume of 175 mm3 at 65 days, GBM50 tumours reached an average volume of 250 mm3 at 61 days (p > 0.05), GBM3 tumours reached an average volume of ~ 150 mm3 at 83 days (p < 0.01), and GBM0.2 tumours reached an average volume of > 200 mm3 at 90 days (p < 0.01) [18].

DNA extraction

Total genomic DNA was extracted from the tumours using the DNeasy Blood & Tissue Kits (Qiagen, CA, USA), according to manufacturer’s protocols with minor modifications. The DNA samples were treated with 3 μL of RNase solution (Qiagen) at room temperature. DNA samples were eluted in 100 μL of autoclaved Milli-Q H2O.

Determination of mtDNA copy number per cell

mtDNA copy number per cell was determined, as previously described [16]. To this extent, concentrations of qPCR products for β-globin and mtDNA were determined against standard curves generated by quantitative real-time PCR (qPCR; Rotor-Gene 3000, Corbett Research, Cambridge, UK) on total DNA purified from the tumours (n = 3 per group; n = 3 replicates per sample). Primer sequences and primer-specific reaction conditions are listed in Additional file 3. mtDNA copy number per cell was calculated using the formula of 2 × NmtDNA/Nβ-globin, where, for NmtDNA and Nβ-globin, N = (qPCR product concentration × 6.023 × 1014)/(qPCR product size in bp × 660).

Immunoprecipitation of methylated DNA sequencing (MeDIP-Seq)

Five microgram of genomic DNA from each of the tumour samples (n = 3 per group) underwent MeDIP, as previously described [70]. Briefly, each DNA sample was sheared into 200–1000 bp using the Covaris Adaptive Focused Acoustics (AFA™) S220 system (Woburn, MA, USA). dsDNA was then denatured to single-stranded DNA, as required for the antibody, by incubation at 95 °C for 10 min. 3 μg of each DNA sample was immunoprecipitated with 2 μg of anti-5mC antibody (Active Motif) with 20 μL per sample of prewashed Dynabeads® Protein G (Thermo Fisher Scientific). The suspension was incubated in 500 μL of IP buffer (100 mM sodium phosphate (pH 7.0); 1.4 M NaCl;0.5% Triton X-100) at 4 °C for 16 h under rotation. The beads were then washed three times with 1 mL of IP buffer and resuspended in 250 μL of proteinase K digestion buffer (50 mM Tris-HCl, pH 8.0; 10 M EDTA, pH 8.0; 1.0% SDS) with 10 μL of proteinase K (20 mg/mL; Bioline, London, UK). The suspension was incubated at 50 °C for 3 h on a thermo-shaker. The supernatant was then collected on a magnetic particle concentrator (Thermo Fisher Scientific). DNA was purified from the supernatant using the QIAquick PCR Purification Kit (Qiagen), according to the manufacturer’s protocol.

MeDIP products from each of the tumour samples underwent library construction using the DNA SMART ChIP-Seq Kit (Clontech, CA, USA), as per the manufacturer’s protocol (No. 021115). Libraries were quality checked using Qubit fluorometric quantitation (Thermo Fisher Scientific), a 2100 Bioanalyzer (Agilent, CA, USA), and qPCR. All libraries were of a similar size (297–335 bp) and quantities (~ 40 ng/μl). A single equimolar pool was made following denaturation of the libraries at 94 °C for 2 min. 12 pM of each library was used for onboard cluster generation and sequencing in 4 lanes of two consecutive Illumina HiSeq 1500 Rapid runs. Cluster densities were within the optimal range of 881–959 k/mm2 (optimal 750–1000 k/mm2). This resulted in around 170 million readable clusters per lane. 95% of the clusters passed filter, and 96% of the reads passed the Illumina sequencing quality score of Q30, which was deemed excellent quality for base-calling. Run quality parameters with the PhiX spike-in had an error rate of 0.12% (expected < 0.5%) and phasing/pre-phasing of < 0.13/< 0.10 (expected < 0.4/< 0.2).

Computational analysis for MeDIP-Seq


Computational analysis was performed on the MeDIP-Seq sequences using the customized analytical pipeline, as described in [31]. Briefly, the sequences were firstly checked for quality and adaptors using the FastQC software (v 0.11.5). The raw sequences were then aligned to the human reference genome GRCh38/hg38 (UCSC) using the Burrows-Wheeler Aligner (BWA) software (version 0.7.16a) in the single-end mode of the BWA-backtrack algorithm with default settings [71]. The output files (*.sam) were then combined and converted to bam files (*.bam) by SAMtools (version 1.4).

Determination of differentially methylated regions (DMRs)

Uniquely mapped reads with mapping scores (MAPQ) ≥ 30 were kept for analytical analysis using the MEDIPS package (version 1.24.0), according to the command codes of the MEDIPS package. Briefly, the reads were cleaned for PCR duplicates that mapped to exactly the same genomic positions by setting the parameter “uniq” to 1. The reads were then extended to 500 nucleotides according to the reference genome without any shift in the genomic locations by setting the parameter “extend” to 500 and “shift” to 0. The coverage of every 100-bp window of the genome was summarized by setting the parameter “window_size” to 100. MeDIP-Seq-specific quality control analysis was performed on each sample, which includes saturation analysis and CpG enrichment analysis. A coupling factor was set up based on the GBM100 sample for normalization. Differentially methylated windows (DMWs) with the total count of reads ≥ 100 (minrowsum = 100) were identified between groups using the ‘edgeR’ method. Statistically significant DMWs were selected for downstream analysis. Continuous significant DMWs were merged as one DMR (Additional file 4).

Annotation of DMRs

DMRs then underwent annotation based on their corresponding genomic regions in the human genome using the ChIPSeeker package (version 3.5) [72]. Super-enhancer annotation was performed using the dbSUPER databases [73]. TAD overlapping analysis was conducted using the TAD datasets available in The ENCODE project [74]. Gene lists were analysed for gene ontology using the Panther Classification System (version 12.0) [75]. The functional classification analysis was used to annotate the functions of genes. The statistical overrepresentation test (Fisher’s exact test with a false discovery rate (FDR) correction) was used to identify the significance of associated biological processes overrepresented amongst the genes.

Region-of-Interest (ROI) analysis

DNA methylation levels at particular regions of the genome, for example CGIs associated with the mtDNA replication factors, were assessed by the built-in ROI analysis in the MEDIPS package. DNA methylation levels of ROIs were expressed as relative methylation scores, as determined by the MEDIPS package. The relative methylation score was developed specifically for MeDIP-Seq in order to normalize methylation scores for regions based on the concept of CpG coupling analysis [76, 77].

Determination of copy number variation (CNV)

Total DNA from the tumour samples (n = 3 per group) was submitted to the Australian Genome Research Facility (AGRF; VIC, AUS) to perform genotyping. The Illumina Human Global Screening Array Beadchip (Illumina) was used, which covers approximately 700 K SNPs throughout the genome. Array data were normalized, clustered and underwent genotype calling using the Genotype Module of GenomeStudio 2.0 (Illumina), according to the user’s manual. The full data report containing log R ratios and B allele frequencies for each probe was then exported and analysed using Nexus 9.0 software (BioDiscovery Inc., CA, USA). The Nexus Copy Number Module was used to identify regions of CN gain and loss, loss of heterozygosity (LOH) for each group and statistically significant CNVs between groups (p value ≤ 0.05). Human genome hg19 was chosen as the reference genome to annotate the genomic locations. Genomic regions were then aligned to human genome hg38, allowing overlapping analysis with the DMRs using the ‘intersect’ function of bedtools (v2.24.0) [78, 79].

Identification of DMR-overlapping differentially expressed genes

RNA-Seq was previously performed on the GBM tumours [17]. The sequencing files were deposited in the NCBI Sequence Read Archive (SRA) under the accession code PRJNA296542 [17]. The sequences were mapped to the human genome (hg19) using the Tophat aligner (v1.3.1) [17, 80]. Differentially expressed genes were determined using the Cufflinks tool (v2.2.1) [17, 80]. Differentially expressed genes (FDR ≤ 0.05) underwent lift-over analysis to match the hg38 human genome assembly using the GALAXY platform ( and then overlapping analysis was performed with the DMRs using the ‘intersect’ function of bedtools (v2.24.0) [78, 79].

Gene expression analysis using real-time quantitative PCR

Total genomic RNA was extracted from the tumour samples (n = 3 per group) using the RNeasy Mini Kit and the QIAshredder (Qiagen, CA, USA), according to manufacturer’s protocols with minor modifications. The RNA samples were treated with DNase I (Qiagen) on column for 20 min. cDNA was synthesized from 1 μg of the total RNA using the Superscript III First-Strand synthesis system (Thermo Fisher Scientific), according to the manufacturer’s instructions.

cDNA products (n = 3 technical replicates per tumour) were assessed using the Rotor-Gene 3000 RT-PCR machine under primer-specific conditions (Additional file 3), as described in [16]. Relative gene expression to the GBM100 tumours was calculated using the ΔΔCT method. OAZ1, 18SrRNA and HPRT1 were used as the housekeeping genes, of which the mean values of expression were used as the internal control for data normalization [31]. Data were represented as the fold change to the GBM100 group (n = 3; mean ± SEM). One-way ANOVA was used to determine statistical significance between the GBM50, GBM3 and GBM0.2 tumours and the GBM100 tumours. Results were plotted using GraphPad Prism 7 (GraphPad Software, Inc., CA, USA).

Gene expression analysis using the Fluidigm platform

Using the same sets of cDNA samples (n = 3 per group; n = 2 replicates per sample), gene expression of targets of interest was assessed using the Fluidigm qPCR array, according to the manufacturer’s instructions. Taqman primers are listed in Additional file 5. Taqman primers were pooled and diluted in C1 DNA suspension buffer to a final concentration for each primer of 180 nM. Each cDNA sample and a non-template control underwent pre-amplification for 14 cycles, according to the manufacturer’s instructions (Quick Reference PN 100-5876 B1). The pre-amplification reaction consisted of the Taqman PreAmp Master mix (Thermo Fisher Scientific) and the pooled Taqman primers. Products were diluted fivefold with C1 DNA suspension buffer. The integrated fluidic circuit controller HX was then used to prime and load the 96.96 Dynamic array plate. 5 μL of each sample was loaded in duplicate into each sample inlet and 5 μL of each Taqman assay (10 x) were loaded into each assay inlet. Real-time qPCR was performed according to the Biomark GE 96.96 Standard v2 protocol. Data were exported using the Fluidigm real-time PCR analysis software (v4.1.1). Relative gene expression was calculated using the same method described above.



mitochondrial DNA


oxidative phosphorylation


polymerase gamma subunit A


polymerase gamma subunit B


mitochondrial helicase


single-stranded DNA binding protein


mitochondrial topoisomerase


mitochondrial transcription factors A


mitochondrial transcription factors B1


mitochondrial transcription factors B2


nuclear respiratory factors 1/2


peroxisome proliferator-activated receptor γ


peroxisome proliferator-activated receptor γ co-activator


Sirtuin family of genes 1-3


oestrogen-related receptors A/B/G


ten-eleven translocation methylcytosine dioxygenases


DNA methyltransferase family


glioblastoma multiforme


methylated DNA immunoprecipitation

GBM100 :

tumours formed from cells possessing 100% of their mtDNA content

GBM50 :

tumours formed from cells possessing 50% of their mtDNA content

GBM3 :

tumours formed from cells possessing 3% of mtDNA content

GBM0.2 :

tumours formed from cells possessing 0.2% of mtDNA content


differentially methylated windows


differentially methylated regions


region of interest


copy number variation






CpG island


topologically associating domains


untranslated regions




Burrows-Wheeler Aligner


false discovery rate


loss of heterozygosity


  1. 1.

    Anderson S, Bankier AT, Barrell BG, de Bruijn MH, Coulson AR, Drouin J, Eperon IC, Nierlich DP, Roe BA, Sanger F, et al. Sequence and organization of the human mitochondrial genome. Nature. 1981;290:457–65.

  2. 2.

    Kasamatsu H, Grossman LI, Robberson DL, Watson R, Vinograd J. The replication and structure of mitochondrial DNA in animal cells. Cold Spring Harb Symp Quant Biol. 1974;38:281–8.

  3. 3.

    Clayton DA. Replication of animal mitochondrial DNA. Cell. 1982;28:693–705.

  4. 4.

    Carrodeguas JA, Theis K, Bogenhagen DF, Kisker C. Crystal structure and deletion analysis show that the accessory subunit of mammalian DNA polymerase gamma, Pol gamma B, functions as a homodimer. Mol Cell. 2001;7:43–54.

  5. 5.

    Takamatsu C, Umeda S, Ohsato T, Ohno T, Abe Y, Fukuoh A, Shinagawa H, Hamasaki N, Kang D. Regulation of mitochondrial D-loops by transcription factor A and single-stranded DNA-binding protein. EMBO Rep. 2002;3:451–6.

  6. 6.

    Korhonen JA, Gaspari M, Falkenberg M. TWINKLE Has 5′ → 3′ DNA helicase activity and is specifically stimulated by mitochondrial single-stranded DNA-binding protein. J Biol Chem. 2003;278:48627–32.

  7. 7.

    Zhang H, Meng LH, Pommier Y. Mitochondrial topoisomerases and alternative splicing of the human TOP1mt gene. Biochimie. 2007;89:474–81.

  8. 8.

    Fisher RP, Clayton DA. Purification and characterization of human mitochondrial transcription factor 1. Mol Cell Biol. 1988;8:3496–509.

  9. 9.

    Falkenberg M, Gaspari M, Rantanen A, Trifunovic A, Larsson NG, Gustafsson CM. Mitochondrial transcription factors B1 and B2 activate transcription of human mtDNA. Nat Genet. 2002;31:289–94.

  10. 10.

    Gleyzer N, Vercauteren K, Scarpulla RC. Control of mitochondrial transcription specificity factors (TFB1M and TFB2M) by nuclear respiratory factors (NRF-1 and NRF-2) and PGC-1 family coactivators. Mol Cell Biol. 2005;25:1354–66.

  11. 11.

    Aquilano K, Vigilanza P, Baldelli S, Pagliei B, Rotilio G, Ciriolo MR. Peroxisome proliferator-activated receptor gamma co-activator 1alpha (PGC-1alpha) and sirtuin 1 (SIRT1) reside in mitochondria: possible direct function in mitochondrial biogenesis. J Biol Chem. 2010;285:21590–9.

  12. 12.

    Chen JQ, Delannoy M, Cooke C, Yager JD. Mitochondrial localization of ERalpha and ERbeta in human MCF7 cells. Am J Physiol Endocrinol Metab. 2004;286:E1011–22.

  13. 13.

    Schreiber SN, Emter R, Hock MB, Knutti D, Cardenas J, Podvinec M, Oakeley EJ, Kralli A. The estrogen-related receptor alpha (ERRalpha) functions in PPARgamma coactivator 1alpha (PGC-1alpha)-induced mitochondrial biogenesis. Proc Natl Acad Sci USA. 2004;101:6472–7.

  14. 14.

    Kong X, Wang R, Xue Y, Liu X, Zhang H, Chen Y, Fang F, Chang Y. Sirtuin 3, a new target of PGC-1alpha, plays an important role in the suppression of ROS and mitochondrial biogenesis. PLoS ONE. 2010;5:e11707.

  15. 15.

    Kelly RD, Mahmud A, McKenzie M, Trounce IA, St John JC. Mitochondrial DNA copy number is regulated in a tissue specific manner by DNA methylation of the nuclear-encoded DNA polymerase gamma A. Nucleic Acids Res. 2012;40:10124–38.

  16. 16.

    Lee W, Johnson J, Gough DJ, Donoghue J, Cagnone GL, Vaghjiani V, Brown KA, Johns TG, St John JC. Mitochondrial DNA copy number is regulated by DNA methylation and demethylation of POLGA in stem and cancer cells and their differentiated progeny. Cell Death Dis. 2015;6:e1664.

  17. 17.

    Lee WT, Cain JE, Cuddihy A, Johnson J, Dickinson A, Yeung KY, Kumar B, Johns TG, Watkins DN, Spencer A, St John JC. Mitochondrial DNA plasticity is an essential inducer of tumorigenesis. Cell Death Discov. 2016;2:16016.

  18. 18.

    Dickinson A, Yeung KY, Donoghue J, Baker MJ, Kelly RD, McKenzie M, Johns TG, St John JC. The regulation of mitochondrial DNA copy number in glioblastoma cells. Cell Death Differ. 2013;20:1644–53.

  19. 19.

    Sun X, St John JC. The role of the mtDNA set point in differentiation, development and tumorigenesis. Biochem J. 2016;473:2955–71.

  20. 20.

    Facucho-Oliveira JM, Alderson J, Splikings EC, Egginton S, St. John JC. Mitochondria DNA replication during differentiation of murine embryonic stem cells. J Cell Sci. 2007;15:4025–34.

  21. 21.

    Facucho-Oliveira JM, St John JC. The relationship between pluripotency and mitochondrial DNA proliferation during early embryo development and embryonic stem cell differentiation. Stem Cell Rev. 2009;5:140–58.

  22. 22.

    El Shourbagy SH, Spikings EC, Freitas M, St John JC. Mitochondria directly influence fertilisation outcome in the pig. Reproduction. 2006;131:233–45.

  23. 23.

    Spikings EC, Alderson J, St John JC. Regulated mitochondrial DNA replication during oocyte maturation is essential for successful porcine embryonic development. Biol Reprod. 2007;76:327–35.

  24. 24.

    Reynier P, May-Panloup P, Chretien MF, Morgan CJ, Jean M, Savagner F, Barriere P, Malthiery Y. Mitochondrial DNA content affects the fertilizability of human oocytes. Mol Hum Reprod. 2001;7:425–9.

  25. 25.

    Santos TA, El Shourbagy S, St John JC. Mitochondrial content reflects oocyte variability and fertilization outcome. Fertil Steril. 2006;85:584–91.

  26. 26.

    Lee WT, St John J. The control of mitochondrial DNA replication during development and tumorigenesis. Ann NY Acad Sci. 2015;1350:95–106.

  27. 27.

    Guo F, Li X, Liang D, Li T, Zhu P, Guo H, Wu X, Wen L, Gu TP, Hu B, et al. Active and passive demethylation of male and female pronuclear DNA in the mammalian zygote. Cell Stem Cell. 2014;15:447–59.

  28. 28.

    Brunner AL, Johnson DS, Kim SW, Valouev A, Reddy TE, Neff NF, Anton E, Medina C, Nguyen L, Chiao E, et al. Distinct DNA methylation patterns characterize differentiated human embryonic stem cells and developing human fetal liver. Genome Res. 2009;19:1044–56.

  29. 29.

    Jones PA, Taylor SM. Cellular differentiation, cytidine analogs and DNA methylation. Cell. 1980;20:85–93.

  30. 30.

    Blaschke K, Ebata KT, Karimi MM, Zepeda-Martinez JA, Goyal P, Mahapatra S, Tam A, Laird DJ, Hirst M, Rao A, et al. Vitamin[thinsp]C induces Tet-dependent DNA demethylation and a blastocyst-like state in ES cells. Nature. 2013;500:222–6.

  31. 31.

    Sun X, Johnson J, St John JC. Global DNA methylation synergistically regulates the nuclear and mitochondrial genomes in glioblastoma cells. Nucleic Acids Res. 2018;46:5977–95.

  32. 32.

    Marie SK, Shinjo SM. Metabolism and brain cancer. Clinics (Sao Paulo). 2011;66(Suppl 1):33–43.

  33. 33.

    Galli R, Binda E, Orfanelli U, Cipelletti B, Gritti A, De Vitis S, Fiocco R, Foroni C, Dimeco F, Vescovi A. Isolation and characterization of tumorigenic, stem-like neural precursors from human glioblastoma. Cancer Res. 2004;64:7011–21.

  34. 34.

    Chu Q, Orr BA, Semenkow S, Bar EE, Eberhart CG. Prolonged inhibition of glioblastoma xenograft initiation and clonogenic growth following in vivo Notch blockade. Clin Cancer Res. 2013;19:3224–33.

  35. 35.

    Spina R, Voss DM, Asnaghi L, Sloan A, Bar EE. Atracurium Besylate and other neuromuscular blocking agents promote astroglial differentiation and deplete glioblastoma stem cells. Oncotarget. 2016;7:459–72.

  36. 36.

    Hess KR, Broglio KR, Bondy ML. Adult glioma incidence trends in the United States, 1977–2000. Cancer. 2004;101:2293–9.

  37. 37.

    Figueroa ME, Abdel-Wahab O, Lu C, Ward PS, Patel J, Shih A, Li Y, Bhagwat N, Vasanthakumar A, Fernandez HF, et al. Leukemic IDH1 and IDH2 mutations result in a hypermethylation phenotype, disrupt TET2 function, and impair hematopoietic differentiation. Cancer Cell. 2010;18:553–67.

  38. 38.

    Turcan S, Rohle D, Goenka A, Walsh LA, Fang F, Yilmaz E, Campos C, Fabius AW, Lu C, Ward PS, et al. IDH1 mutation is sufficient to establish the glioma hypermethylator phenotype. Nature. 2012;483:479–83.

  39. 39.

    Yan H, Parsons DW, Jin G, McLendon R, Rasheed BA, Yuan W, Kos I, Batinic-Haberle I, Jones S, Riggins GJ, et al. IDH1 and IDH2 mutations in gliomas. N Engl J Med. 2009;360:765–73.

  40. 40.

    Cohen AL, Holmen SL, Colman H. IDH1 and IDH2 mutations in gliomas. Curr Neurol Neurosci Rep. 2013;13:345.

  41. 41.

    Robinson MD, Stirzaker C, Statham AL, Coolen MW, Song JZ, Nair SS, Strbenac D, Speed TP, Clark SJ. Evaluation of affinity-based genome-wide DNA methylation data: effects of CpG density, amplification bias, and copy number variation. Genome Res. 2010;20:1719–29.

  42. 42.

    Nothjunge S, Nuhrenberg TG, Gruning BA, Doppler SA, Preissl S, Schwaderer M, Rommel C, Krane M, Hein L, Gilsbach R. DNA methylation signatures follow preformed chromatin compartments in cardiac myocytes. Nat Commun. 2017;8:1667.

  43. 43.

    Haseeb A, Makki MS, Haqqi TM. Modulation of ten-eleven translocation 1 (TET1), Isocitrate Dehydrogenase (IDH) expression, alpha-Ketoglutarate (alpha-KG), and DNA hydroxymethylation levels by interleukin-1beta in primary human chondrocytes. J Biol Chem. 2014;289:6877–85.

  44. 44.

    Calvert AE, Chalastanis A, Wu Y, Hurley LA, Kouri FM, Bi Y, Kachman M, May JL, Bartom E, Hua Y, et al. Cancer-associated IDH1 promotes growth and resistance to targeted therapies in the absence of mutation. Cell Rep. 2017;19:1858–73.

  45. 45.

    Tan AS, Baty JW, Dong LF, Bezawork-Geleta A, Endaya B, Goodwin J, Bajzikova M, Kovarova J, Peterka M, Yan B, et al. Mitochondrial genome acquisition restores respiratory function and tumorigenic potential of cancer cells without mitochondrial DNA. Cell Metab. 2015;21:81–94.

  46. 46.

    Dong LF, Kovarova J, Bajzikova M, Bezawork-Geleta A, Svec D, Endaya B, Sachaphibulkij K, Coelho AR, Sebkova N, Ruzickova A, et al. Horizontal transfer of whole mitochondria restores tumorigenic potential in mitochondrial DNA-deficient cancer cells. Elife. 2017;6:e22187.

  47. 47.

    Yeung KY, Dickinson A, Donoghue JF, Polekhina G, White SJ, Grammatopoulos DK, McKenzie M, Johns TG, St John JC. The identification of mitochondrial DNA variants in glioblastoma multiforme. Acta Neuropathol Commun. 2014;2:1.

  48. 48.

    Heller S, Schubert S, Krehan M, Schafer I, Seibel M, Latorre D, Villani G, Seibel P. Efficient repopulation of genetically derived rho zero cells with exogenous mitochondria. PLoS ONE. 2013;8:e73207.

  49. 49.

    Deaton AM, Bird A. CpG islands and the regulation of transcription. Genes Dev. 2011;25:1010–22.

  50. 50.

    Yang X, Han H, De Carvalho DD, Lay FD, Jones PA, Liang G. Gene body methylation can alter gene expression and is a therapeutic target in cancer. Cancer Cell. 2014;26:577–90.

  51. 51.

    Bell RE, Golan T, Sheinboim D, Malcov H, Amar D, Salamon A, Liron T, Gelfman S, Gabet Y, Shamir R, Levy C. Enhancer methylation dynamics contribute to cancer plasticity and patient mortality. Genome Res. 2016;26:601–11.

  52. 52.

    Shen L, Zhang Y. Enzymatic analysis of Tet proteins: key enzymes in the metabolism of DNA methylation. Methods Enzymol. 2012;512:93–105.

  53. 53.

    Scourzic L, Mouly E, Bernard OA. TET proteins and the control of cytosine demethylation in cancer. Genome Med. 2015;7:9.

  54. 54.

    Cuyas E, Fernandez-Arroyo S, Verdura S, Garcia RA, Stursa J, Werner L, Blanco-Gonzalez E, Montes-Bayon M, Joven J, Viollet B, et al. Metformin regulates global DNA methylation via mitochondrial one-carbon metabolism. Oncogene. 2018;37:963–70.

  55. 55.

    Chen H, Wang J, Liu Z, Yang H, Zhu Y, Zhao M, Liu Y, Yan M. Mitochondrial DNA depletion causes decreased ROS production and resistance to apoptosis. Int J Mol Med. 2016;38:1039–46.

  56. 56.

    Kelly RD, Sumer H, McKenzie M, Facucho-Oliveira J, Trounce IA, Verma PJ, St John JC. The effects of nuclear reprogramming on mitochondrial DNA replication. Stem Cell Rev. 2013;9:1–15.

  57. 57.

    Hance N, Ekstrand MI, Trifunovic A. Mitochondrial DNA polymerase gamma is essential for mammalian embryogenesis. Hum Mol Genet. 2005;14:1775–83.

  58. 58.

    Sajadian SO, Tripura C, Samani FS, Ruoss M, Dooley S, Baharvand H, Nussler AK. Vitamin C enhances epigenetic modifications induced by 5-azacytidine and cell cycle arrest in the hepatocellular carcinoma cell lines HLE and Huh7. Clin Epigenetics. 2016;8:46.

  59. 59.

    Liu M, Ohtani H, Zhou W, Orskov AD, Charlet J, Zhang YW, Shen H, Baylin SB, Liang G, Gronbaek K, Jones PA. Vitamin C increases viral mimicry induced by 5-aza-2′-deoxycytidine. Proc Natl Acad Sci USA. 2016;113:10238–44.

  60. 60.

    Berezovsky AD, Poisson LM, Cherba D, Webb CP, Transou AD, Lemke NW, Hong X, Hasselbach LA, Irtenkauf SM, Mikkelsen T, deCarvalho AC. Sox2 promotes malignancy in glioblastoma by regulating plasticity and astrocytic differentiation. Neoplasia. 2014;16:193–206, 206 e119–125.

  61. 61.

    Muona M, Berkovic SF, Dibbens LM, Oliver KL, Maljevic S, Bayly MA, Joensuu T, Canafoglia L, Franceschetti S, Michelucci R, et al. A recurrent de novo mutation in KCNC1 causes progressive myoclonus epilepsy. Nat Genet. 2015;47:39–46.

  62. 62.

    Aref-Eshghi E, Rodenhiser DI, Schenkel LC, Lin H, Skinner C, Ainsworth P, Pare G, Hood RL, Bulman DE, Kernohan KD, et al. Genomic DNA methylation signatures enable concurrent diagnosis and clinical genetic variant classification in neurodevelopmental syndromes. Am J Hum Genet. 2018;102:156–74.

  63. 63.

    Hoischen A, Krumm N, Eichler EE. Prioritization of neurodevelopmental disease genes by discovery of new mutations. Nat Neurosci. 2014;17:764–72.

  64. 64.

    Reis K, Fransson A, Aspenstrom P. The Miro GTPases: at the heart of the mitochondrial transport machinery. FEBS Lett. 2009;583:1391–8.

  65. 65.

    Kay L, Pienaar IS, Cooray R, Black G, Soundararajan M. Understanding miro GTPases: implications in the treatment of neurodegenerative disorders. Mol Neurobiol. 2018;55:7352–65.

  66. 66.

    Lottaz C, Beier D, Meyer K, Kumar P, Hermann A, Schwarz J, Junker M, Oefner PJ, Bogdahn U, Wischhusen J, et al. Transcriptional profiles of CD133 + and CD133- glioblastoma-derived cancer stem cell lines suggest different cells of origin. Cancer Res. 2010;70:2030–40.

  67. 67.

    Campbell AM, Zagon IS, McLaughlin PJ. Astrocyte proliferation is regulated by the OGF-OGFr axis in vitro and in experimental autoimmune encephalomyelitis. Brain Res Bull. 2013;90:43–51.

  68. 68.

    Lee WT, Sun X, Tsai TS, Johnson JL, Gould JA, Garama DJ, Gough DJ, McKenzie M, Trounce IA, St John JC. Mitochondrial DNA haplotypes induce differential patterns of DNA methylation that result in differential chromosomal gene expression patterns. Cell Death Discov. 2017;3:17062.

  69. 69.

    Martin JL, Brown CE, Matthews-Davis N, Reardon JE. Effects of antiviral nucleoside analogs on human DNA polymerases and mitochondrial DNA synthesis. Antimicrob Agents Chemother. 1994;38:2743–9.

  70. 70.

    Weber M, Davies JJ, Wittig D, Oakeley EJ, Haase M, Lam WL, Schubeler D. Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. Nat Genet. 2005;37:853–62.

  71. 71.

    Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–60.

  72. 72.

    Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015;31:2382–3.

  73. 73.

    Khan A, Zhang X. dbSUPER: a database of super-enhancers in mouse and human genome. Nucleic Acids Res. 2016;44:D164–71.

  74. 74.

    ENCODE Project Consortium. The ENCODE (ENCyclopedia Of DNA Elements) Project. Science. 2004;306:636–40.

  75. 75.

    Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, Thomas PD. PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2017;45:D183–9.

  76. 76.

    Chavez L, Jozefczuk J, Grimm C, Dietrich J, Timmermann B, Lehrach H, Herwig R, Adjaye J. Computational analysis of genome-wide DNA methylation during the differentiation of human embryonic stem cells along the endodermal lineage. Genome Res. 2010;20:1441–50.

  77. 77.

    Lienhard M, Grimm C, Morkel M, Herwig R, Chavez L. MEDIPS: genome-wide differential coverage analysis of sequencing data derived from DNA enrichment experiments. Bioinformatics. 2014;30:284–6.

  78. 78.

    Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

  79. 79.

    Quinlan AR. BEDTools: the Swiss-army tool for genome feature analysis. Curr Protoc Bioinform. 2014;47:11.12.1–34.

  80. 80.

    Ghosh S, Chan CK. Analysis of RNA-seq data using TopHat and cufflinks. Methods Mol Biol. 2016;1374:339–61.

Download references

Authors’ contributions

XS designed and performed the experiments, analysed the data and drafted the manuscript. JCSJ conceived the work, designed and coordinated the experiments, analysed the data, drafted the manuscript and obtained funding for the work. Both authors read and approved the final manuscript.


We are grateful to Dr. Trevor Wilson and Dr. Selva Kumari Ramasubramanian, The Medical Genomics Facility, Monash Health Translation Precinct, for assistance in performing MeDIP-Seq and the Fluidigm real-time qPCR array, respectively. We are also grateful to BioDiscovery, Inc for the use of Nexus Copy Number software in our data analysis.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

The MeDIP-Seq datasets supporting the conclusions of this article are available in Sequence Read Archive (SRA) and are accessible through accession number PRJNA397951. The CNV datasets supporting the conclusions of this article are available in the Gene Expression Omnibus (GEO) and are available under accession number GSE98366. The RNA-Seq datasets supporting the conclusions of this article are available in SRA and are obtainable under accession number PRJNA296542, as detailed in [17].

Consent for publication

Not applicable.

Ethics approval and consent to participate

The animal work was approved by the Animal Ethics Committee, Monash University, Approval Number: MMCA/2011/76.


This work was supported by Hudson Institute of Medical Research Discretionary Funds, and the Victorian Government’s Operational Infrastructure Support Program. XS was supported by an Australian Postgraduate Award.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Author information

Correspondence to Justin C. St John.

Additional files


Additional file 1. MeDIP-Seq specific QC results determined by the MEDIPS package.


Additional file 2. Summary of CNV regions identified in each cohort of tumours and overlapping with DMRs using the Nexus 9.0 software.


Additional file 3. Primer pairs for real time PCR.


Additional file 4. DMRs identified using MeDIP-Seq and MEDIPS package.


Additional file 5. Taqman assays used in the Fluidigm qPCR arrays.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Sun, X., St John, J.C. Modulation of mitochondrial DNA copy number in a model of glioblastoma induces changes to DNA methylation and gene expression of the nuclear genome in tumours. Epigenetics & Chromatin 11, 53 (2018).

Download citation


  • Mitochondrial DNA
  • DNA methylation
  • Gene expression
  • Tumorigenesis
  • POLG