Human ovarian cancer cells SKOV3 were a gift from Kristian Helin’s lab (Copenhagen). Human embryonic kidney cells HEK293 and SKOV3 cells were maintained in DMEM with 20% FBS and 1% penicillin-streptomycin in an atmosphere of 5% CO₂ at 37°C.

Lonza human skeletal muscle myoblasts (Lonza) were maintained in DMEM/F-12, GlutaMAX^TM (Gibco, Life Technologies) supplemented with 20% Fetal Bovine Serum (FBS, Sigma-Aldrich) and 1% penicillin-streptomycin (Invitrogen) in 5% CO₂ at 37°C. At 70–80% confluency, differentiation was induced by reducing the FBS concentration to 2%. Experiments were conducted when cells were fully differentiated myotubes, after 5–7 days in differentiation medium.

Wild type C57BL/6J male mice (Taconic Biosciences, Denmark) were 6–8 weeks of age at delivery and single-housed for acclimatization 1 week prior to the start of the experiment. All experiments were conducted in accordance with national Danish guidelines (amendment #1306 of November 23, 2007) as approved by the Danish Animal Experiments Inspectorate (#2014-15-2934-01,027). Mouse muscle tissues were collected from donor mice of the previously published experiment (Jensen et al., 2016).

Mitochondria from human myotubes were isolated with an adapted protocol described elsewhere (Mercer et al., 2011). Cells were trypsinized with Trypsin-2.5% EDTA and resuspended in PBS. Cells were resuspended in swelling buffer (10 mM NaCl, 1.5 mM MgCl₂ and 10 mM Tris-HCl, pH 7.5 all from Sigma-Aldrich) and were allowed to swell for 5 min on ice. Cells were then homogenized seven times with a 7-ml glass homogenizer and sucrose concentration was adjusted to 250 mM. Nuclei were sedimented by centrifugation at 1300 g for 3 min at 4°C. The supernatant containing the mitochondria was transferred to another tube and centrifuged at 15,000 g for 10 min at 4°C. The mitochondrial pellet was centrifuged for an additional 15 min at 15,000 g at 4°C. The mitochondrial pellet was resuspended in PBS and mtDNA was extracted using DNeasy Blood and Tissue Kit (Qiagen) and quantified using Qubit dsDNA HS Kit (Life Technologies).

Mitochondria from mouse gastrocnemius were isolated using an adapted protocol described elsewhere (Lanza and Nair, 2009). The muscle was placed immediately into ice-cold muscle homogenization buffer (100 mM KCl, 50 mM Tris-HCl, 5 mM MgCl₂, 1.8 mM ATP, and 1 mM EDTA pH 7.2, all from Sigma-Aldrich). Fat and connective tissues were removed, and the muscle was minced. Muscle pieces incubated for 2 min in 1 mL of protease medium [1 mL homogenization buffer and 5.66 mg protease from Bacillus licheniformis, 10.6 U/mg (Sigma-Aldrich)], washed twice with 3 mL of homogenization buffer, and transferred to a homogenizer containing 4 mL of homogenization buffer. The muscle was homogenized using a motor-driven homogenizer for 5 min at 150 rpm with slow vertical plunger movements. The homogenate was then centrifuged at 900 g for 5 min at 4°C. The supernatant containing the mitochondria was collected, and centrifuged a first time at 900 g for 5 min 4°C and a second time at 10,000 g for 5 min at 4°C to pellet the mitochondria. The pellet was washed once with homogenization buffer, and centrifuged for 5 min at 9,000 g at 4°C. The mitochondrial pellet was finally resuspended in PBS and mtDNA was extracted using DNeasy Blood and Tissue Kit (Qiagen). RNase A (Qiagen) treatment was included in the protocol.

Genomic DNA was either untreated or treated with BamHI restriction enzyme (New England Biolabs, cat # R0136) to obtain circular or linearized mtDNA. BamHI recognize and cuts G/GATCC, which is found only at one site in the human mitochondrial DNA at position 14258. The DNA was treated for 4 h at 37°C according to the recommendations from the manufacturer.

Genomic DNA was bisulfite converted using EZ DNA methylation-lightning kit (Zymo Research, United States) according to the manufacturer’s instructions. The converted DNA was amplified using bisulfite converted primers designed using the online primer design programs BiSearch (Tusnády et al., 2005; Arányi et al., 2006) or Methprimer (Li and Dahiya, 2002) (Table 1). Primers for the D-loop, tRNA-f+12S and ND5 targeted the heavy strand, whereas primers for 16S and CYTB targeted the light strand. Primers were checked for specificity using BiSearch to ensure that only mtDNA was amplified. mtDNA amplification was performed using the HotStar Taq plus DNA polymerase kit (Qiagen) according to the manufacturer’s instructions. Briefly, 100 ng converted DNA, 500 μM sense and anti-sense primers and HotStar Taq polymerase 7.5 units/reaction were mixed in a total volume of 50 μl. The PCR was run with the following conditions: 5 min 95°C (60 s 94°C; 60 s 53–55°C; 60 s 72°C) × 35 cycles; 10 min 72°C. Primers were either used separately or multiplexed. When multiplexing primers 200 ng converted DNA was used and the following cycling conditions: 5 min 95°C (60 s 94°C; 90 s 53–55°C; 90 s 72°C) × 35 cycles; 10 min 72°C. PCR products were subjected to gel electrophoresis on 2% agarose (Lonza) and extracted using QIAquick Gel Extraction Kit (Qiagen).

Library preparation of bisulfite-converted PCR products was performed using NEBNext Ultra DNA Library Kit for Illumina (New England Biolabs) according to the manufacturer’s instructions. In short, 100 ng of PCR product was end-repaired and ligated to adaptors. During PCR amplification indexes from NEBNext Singleplex Kit (New England Biolabs) were introduced and libraries were purified using AMPure beads (Beckman Coulter). Libraries were subjected to 150 bp paired-end sequencing using an Illumina MiSeq instrument (Illumina). Total number of reads was approximately 15 million.

Mitochondria were isolated and mtDNA was purified as described above. Whole mitochondrial genome bisulfite sequencing (WMGBS) was performed using EpiGnome^TM Methyl-Seq Kit (Epicentre, Illumina) according to the manufacturer’s protocol. Mitochondrial DNA was sonicated 2 × 10 min 30 Hz using Bioruptor (Diagenode) and bisulfite converted using EZ DNA Methylation-Lightning kit (Zymo Research). 50 ng converted mtDNA was used for library preparation. Random hexamers were added to converted mtDNA following first strand synthesis. After the addition of tag sequences to the 3′ and 5′ ends, a PCR was performed to amplify libraries for 2nd strand synthesis and to add both adaptor and indexes. Libraries were purified using AMPure beads (Beckman coulter) and quality-controlled on a Bioanalyzer 2100 (Agilent Technologies). Libraries were subjected to 150 bp paired-end sequencing using Illumina MiSeq instrument (Illumina).

Whole mitochondrial genome reads and targeted regions were processed similarly: reads were preprocessed with Trim Galore v0.4.1 & Cutadapt v1.12 (Martin, 2011) with the --paired and --trim1 flags set. Due to biases introduced by the primers (visible in an M-bias plot, data not shown) the first nucleotides of each read were trimmed, whole genome reads had 10 nucleotides removed while reads from targeted regions had 2 removed. Preprocessed reads were aligned to the hg38 or mm10 genome depending on dataset using Bismark v0.16.3 (Krueger and Andrews, 2011) assisted by Bowtie2 v2.2.8 (Langmead and Salzberg, 2012) and with the --dovetail flag set. Whole mitochondrial genome reads were de-duplicated using the Bismark de-duplication module. Reads mapping to the mitochondrial genome were extracted using samtools v1.3.1 (Li et al., 2009). Finally, methylation individual methylation levels were extracted using Bismark methylation with the flags --cytosine_report and --CX_context.

Further analysis and visualization was done in the R statistical environment (R Development Core Team, 2014), using the packages data.table (Anon, 2017) and ggplot2 (Wickham, 2009). Testing for the effect of linearizing mtDNA was done using a sign-test and adjusted by Bonferroni correction. Smoothed conditional means were calculated using the “gam” method and a span of 0.2. All correlations calculated are Pearson correlations.

Statistical analysis was performed using R statistical environment software. In whole genome bisulfite sequencing correlations between estimated methylation and sequencing depth were based on the ∼6,000 to ∼12,000 mitochondrial cytosines covered by the experiment with an estimated methylation greater than 0. The Pearson correlation coefficient was calculated between the logarithm of the number of reads covering a cytosine and the logarithm of the estimated methylation. A P-value ≤ 0.05 was considered statistically significant.

To map methylation patterns across the whole mitochondrial genome at single base-pair resolution, we performed whole mitochondrial genome bisulfite sequencing (WMGBS). Methylation of mtDNA has been detected in mouse skeletal muscle tissue (Wong et al., 2013) but not in human Epithelial Kidney cells (HEK293) (Hong et al., 2013). Thus, we chose mouse muscle tissue, human muscle and HEK293 cells to study mtDNA methylation. While cytosine unconversion rate in WMGBS was detected up to 100% at certain positions, suggesting 100% methylation at a few specific mtDNA sites, it was striking to observe that cytosine unconversion rate was negatively associated with sequencing depth: in human cell lines and in mouse skeletal muscle, mtDNA methylation at specific cytosine residues was higher when these residues had a low sequencing coverage (Figure 1 and Supplementary Figures 1A–E). Indeed, calculation of the Pearson correlation coefficient between mtDNA methylation and sequencing depth returned a very tight negative correlation (R = -0.907, P < 2.2E-16; Figure 2). The fact that mtDNA from HEK293 cells had an overall higher sequencing depth than other cells or tissues and was associated to a lower detection of methylation levels, further supports that sequencing depth can influence the estimation of mtDNA methylation (Figure 1).

Given the inevitable contamination of mtDNA preparation with nuclear DNA (nDNA), nuclear insertions of the mitochondrial genome, named NUclear MiTochondrial segments (NUMTs) (Ramos et al., 2011), represent a potential source of analysis bias after WMGBS. While short NUMTs and NUMTs with low sequence similarity are easy to computationally distinguish from mtDNA, longer NUMTs or NUMTs with sequence similarities could influence the quantification of mtDNA methylation. However, we found no correlation between methylation levels and the length of NUMTs (R = 0.06, P = 0.11; Supplementary Figure 2A), and only a marginal correlation between methylation and NUMT alignment score, a score established to represent how well each NUMTs align to mtDNA (Ramos et al., 2011) (R = -0.1, P = 0.006; Supplementary Figure 2B). Bisulfite conversion reduces the sequence complexity of unmethylated DNA and thus, distinguishing NUMTs from mtDNA could be more difficult after bisulfite conversion. To ensure that we could distinguish NUMTs in our analysis, we extracted all possible 300 bp reads from both methylated and unmethylated NUMTs and aligned them back onto the genome using the same setting as in the main analysis. Of the 900.000+ possible reads, none of them aligned to the mitochondrial genome. These results indicate that in our mtDNA preparations, contamination of nDNA was not at the origin of the methylation signal.

The difference in sequencing depth across the mtDNA that we detected in our WMGBS analysis could originate from the secondary structure of mtDNA. The mtDNA genome can be organized in coiled and supercoiled secondary structures (Kolesar et al., 2013). During sonication of mtDNA, a step of the WGBS protocol, certain mtDNA fragments that are less tightly engaged in the mtDNA supercoiled structure have, in theory, the potential to be preferably freed and ligated to sequencing adapters. This would result in an overrepresentation of mtDNA fragments that do not engage in supercoiling. We tested the hypothesis that mtDNA supercoils can be a source of overestimation of methylation by subjecting DNA to digestion with the BamHI restriction enzyme prior to targeted bisulfite sequencing to remove any secondary structure and linearize mtDNA. We focused our analysis on five mitochondrial regions near the origin of replication and transcription start site for the heavy strand: the displacement Loop (D-Loop), tRNA phenylalanine and 12S ribosome (tRNA-F+12S), 16S ribosome (16S), NADH dehydrogenase 5 (ND5) and cytochrome b (CYTB) mRNA-encoding gene (sequencing primers are shown in Table 1) in two human cell lines; In human myotubes and the ovarian cancer cell line SKOV3. mtDNA methylation has previously been detected in SKOV3 cells (van der Wijst et al., 2017). In the D-loop (position 6–298), we found methylation levels from 0–4.8% in Lonza and 0–16.9% in SKOV3 cells with undigested mtDNA. In contrast, the same region after BamHI digestion showed substantial lower methylation levels (0–0.6%, 0–1.1% P = 2.02E-41, P = 1.40E-24). Similarly, in the position 279–458 of the D-loop, we detected 5.5–15.1% and 0.5–1% methylation in undigested vs. digested mtDNA from Lonza, respectively, (P = 5.00E-12) whereas for SKOV3 cells, we detected 19–32.3% and 0.8–1.5% (P = 8.20E-08). In the regions encoding tRNA phenylalanine and the 12S ribosomal RNA and 16S ribosomal RNA, the methylation percentage dropped from 0–7.9% to 0–1.1% (P = 9.22E-15) and 0–0.9% to 0–0.6% (P = 6.50E-05) in Lonza and 0–11.1% to 0–1.1% in SKOV3 cells (P = 1.30E-9), respectively. Finally, in the ND5 and CYTB regions the drop was from 0–3.8% to 0–1.3% (P = 1.7E-05) and 0–1.2% to 0–0.5% (P = 2.95E-09) after digestion for LONZA and from 0–7.9% to 0–1.1% (P = 1.0E-12) and 0–1.8% to 0–0.8% (P = 6.68E-13) for SKOV3 cells, respectively (Figure 3). These results suggest that mtDNA secondary structure impairs bisulfite conversion and infer that the secondary structure of mtDNA is a source of overestimation of methylation levels in bisulfite sequencing.