CpG methylation patterns of human mitochondrial DNA

The epigenetic modification of mitochondrial DNA (mtDNA) is still in controversy. To clarify this point, we applied the gold standard method for DNA methylation, bisulfite pyrosequencing, to examine human mtDNA methylation status. Before bisulfite conversion, BamHI was used to digest DNA to open the loop of mtDNA. The results demonstrated that the linear mtDNA had significantly higher bisulfite conversion efficiency compared with circular mtDNA. Furthermore, the methylation values obtained from linear mtDNA were significantly lower than that of circular mtDNA, which was verified by SEQUENOM MassARRAY. The above impacts of circular structure were also observed in lung DNA samples but not in saliva DNA samples. Mitochondrial genome methylation of blood samples and saliva samples from 14 unrelated individuals was detected. The detected regions covered 83 CpG sites across mtDNA including D-loop, 12 S rRNA, 16 S rRNA, ND1, COXI, ND3, ND4, ND5, CYTB. We found that the average methylation levels of nine regions were all less than 2% for both sample types. In conclusion, our findings firstly show that the circular structure of mtDNA affects bisulfite conversion efficiency, which leads to overestimation of mtDNA methylation values. CpG methylation in human mtDNA is a very rare event at most DNA regions.

Nowadays, several innovative approaches for measuring DNA methylation have been established based on different techniques. For instance, liquid chromatography-electrospray ionization tandem mass spectrometry (LC-ESI-MS/MS) offers a fundamental tool for global quantitation of DNA methylation 10 . By applying this method, human global mtDNA was detected to undergo methylation 11 . ELISA was also used to verify the existence of 5 mC in mtDNA 12,13 . Gene specific methylation approaches are mainly based on affinity enrichment, restriction enzymes digestion or bisulfite conversion 10 . Using affinity-based methods such as 5 mC and 5 hmC immunoprecipitation, both 5 mC and 5 hmC modification were found to be present in D-loop region of mtDNA in human and mouse blood and cultured cells 14,15 . Although affinity-based methods are fast and efficient, they cannot provide the methylation status of individual CpG dinucleotides. Bisulfite conversion coupled with sequencing is the best way to evaluate the methylation status with single base resolution. Recently, several studies have evaluated the methylation of mtDNA by using bisulfite sequencing or pyrosequencing or specific PCR, most of the studies reported that mtDNA methylation could be detected within a level of 2% to 18% in D-loop region or 12 s rRNA or 16 s rRNA in human peripheral blood, cord blood, placenta, and mouse brain, liver, and testes [16][17][18][19] . On the contrary, two studies showed that methylation of mtDNA was a rare event or absent in the detected specific regions in human cell lines and primary cells 20,21 .
These discordant results let us to discover the real reasons for such a dispute and to map human mitochondrial genome methylation. Although a recent study by Ghosh et al described a comprehensive map of methyl cytosines across human mitochondrial genome by analyzing 39 MeDIP datasets and a bisulfite sequencing dataset, it did not provide any information on the sites or frequency of CpG methylation of the targets 22 . Bisulfite pyrosequencing technology is widely recognized as the gold standard of DNA methylation analysis, which can get precise quantification of each CpG site, and has higher accuracy through internal quality control monitoring bisulfite conversion efficiency 23,24 . Therefore, this study applied bisulfite pyrosequencing technology to map the methylation profile of human mtDNA.

Results and Discussion
Circular structure of mtDNA affects bisulfite conversion. To observe whether the ring structure of mtDNA affects the bisulfite conversion efficiency and pyrosequencing results, genomic DNA was treated with BamHI before bisulfite conversion. Four sequences MT1, MT3, MT10, MT13 located in D-loop, 12 s rRNA, ND1, and COX I (Table 1) were chosen to analyze the loop structure effect on bisulfite conversion. Through the detection of 154 blood samples, we found the sequencing success rate of four sequences MT1, MT3, MT10, MT13 in linear mtDNA group were 78.07%, 97.32%, 97.73%, 90.26%, while in circular mtDNA group were 80.95%, 84.62%, 64.52%, 47.89%, respectively. According to statistical analysis, except MT1 (P = 0.621), the success rate of the remaining three sequences between linear mtDNA group and circular mtDNA group had significant differences, and the success rate of circular mtDNA group were significantly lower than that of linear mtDNA group (P = 0.002 for MT3, P < 0.001for MT10 and MT13, Fig. 1).
The main reasons leading to failure in bisulfite pyrosequencing include incomplete bisulfite conversion that can be evaluated by internal quality control during pyrosequencing, low signal-to-noise ratio or signal loss caused by insufficient PCR products, and uncertain reference sequence pattern due to stochastic effect. Through further analyzing the causes for sequencing failure, we found that in the circular mtDNA group, the proportion of the unpassed inner quality control in all of the four sequences except MT3, were more than 80%. While in the linear mtDNA group, the main reason for sequencing failure was low signal-to-noise ratio or signal loss (Table 2), which were usually related to the experiment operation or caused by some stochastic effect. Failure of inner quality control was shown as higher C/T ratio than threshold in non-CpG cytosine, indicating that bisulfite conversion was not complete. Therefore, the above results showed that the circular structure of mtDNA in blood samples obviously affected the bisulfite conversion process, and the conversion efficiency was significantly decreased compared with that of linear mtDNA treated by BamHI. Methylation modification in mtDNA is still in dispute, although numerous innovative approaches for measuring DNA methylation have been used to measure it. This study showed that the circular structure of mtDNA was a possible reason for the discordant results obtained by using bisulfite-based methods. To our knowledge, it is the first time to compare the methylation analysis differences between primary mtDNA with loop structure and linear mtDNA treated with restriction endonuclease. We found that the circular structure of mtDNA hampered the bisulfite conversion process to some extent. This is an important finding that has previously been neglected.
Overestimation of CpG methylation values due to uncompleted bisulfite conversion. Although the sequencing success rate of the circular mtDNA group in blood samples was greatly lower, more than 47% of the samples still succeeded in pyrosequencing. Therefore, we compared the methylation quantification values between the two groups. The results showed that the average methylation levels of MT1, MT3, MT10, MT13 were 11.65% ± 0.29%, 3.58% ± 0.27%, 5.68% ± 0.20%, 7.21% ± 0.27% in circular mtDNA group, while they were 1.16% ± 0.09%, 1.68% ± 0.10%, 0.32% ± 0.05%, 0.78% ± 0.09%, respectively in linear mtDNA group. Wilcoxon signed-rank test showed that significant differences exist between the two groups (P < 0.001 for all of the CpG sites), and the methylation degree in linear mtDNA group was significantly lower than that of circular mtDNA group (Fig. 2). For the same samples tested by both methods and for which methylation values were obtained, paired statistic test was performed and the results demonstrated that the D values within two methods of all the 17 CpG sites except CpG39 were statistically significant (Table 3), which further confirmed that the methylation values obtained from the circular mtDNA were significantly higher than that of the linear ones.
Therefore, if we apply bisulfite-based methods to analyze mtDNA methylation without opening loop before conversion, the results are unreliable. Our findings remind us that even pyrosequencing is performed without opening-loop treatment, the methylation values tend to be overestimated. Bisulfite-based sequencing has been used to examine cytosine modification for 20 years and remains the "gold standard" for DNA methylation analysis 25 . After bisulfite treatment, unmodified cytosine can be converted to uracil, but modified 5 mC and 5 hmC cannot be converted. Methylation modification is then changed to T/C variation, which can be detected and the methylation levels can be quantified by calculating the C/T ratio. So, if this process is not efficient, some unmodified cytosine will be not converted and remain as cytosine. As a result, C/T ratio will be increased and the methylation values will be overestimated. Some studies using bisulfite pyrosequencing to measure the methylation level of mtDNA obtained higher methylation values than ours. For example, Byun et al. 16 measured methylation by bisulfite pyrosequencing in three mtDNA regions in human peripheral blood samples. They reported that the mean methylation level was 5.06% for 12 s rRNA and 2.38% for D-loop region, while ours are 1.68% for 12 s rRNA and 1.16% for D-loop region in the same kind of sample. Besides, detected by bisulfite-based methods, mean methylation degree has been reported to be 4.0% in D-loop region in human cord blood 18 , 3.7% in human placenta 18 , and 2 ~ 9% in mouse brain, liver and testes 17 . And for 12 s rRNA, the methylation value has been reported to be 11.7% in human cord blood 18 , 9.5% in placenta 18 . These methylation values might be overestimated due to incomplete bisulfite conversion.
We noticed that the study conducted by Hong et al. 20 digested DNA with HindIII before bisulfite conversion. HindIII specifically recognizes AAGCTT and cuts between A/A. Three HindIII restriction sites exist in mtDNA. Therefore mtDNA will be cut into linear DNA after HindIII digestion. They then used bisulfite sequencing to detect methylation level of 12 Sr RNA, 16 S rRNA, COII and ATP6 in mitochondrial DNA in human blood cells and HCT116 cell line. The results showed that the methylation degree of these four regions Figure 1. The circular structure of mtDNA affects the pyrosequencing pass ratio in blood samples but not in saliva samples. In blood samples, pyrosequencing pass ratio in linear mtDNA group was significantly higher than that of circular mtDNA group. In saliva samples, there were no significant differences between two groups. * P < 0.01 vs linear mtDNA group.
were no more than 0.66%, and confirmed the absence of CpG methylation in human mtDNA, which was consistent with our studies. Therefore, we can conclude clearly that the loop structure of mtDNA must be opened before bisulfite treatment if using any bisulfite-based methods to measure mtDNA methylation, at least for blood samples, or the methylation values will be overestimated. It has been reported that bisulfite treatment might result in the fragmentation of DNA. The conditions necessary for complete conversion, such as long incubation times, elevated temperature, and high bisulfite concentration, can lead to the degradation of the  incubated DNA 26 . However, since mtDNA is more stable and more anti-degraded than nuclear DNA due to the circular structure, we speculate that the intactness of mtDNA conformation might be reserved more or less. It should be confirmed by further study.

Validation by SEQUENOM EpiTYPER MassARRAY.
To validate the results obtained by pyrosequencing, we reanalyzed the mtDNA methylation status using SEQUENOM EpiTYPER MassARRAY. SEQUENOM EpiTYPER MassARRAY assay is one of the first truly high-throughput and quantitative methods to analyze DNA methylation. EpiTYPER was first introduced by Ehrich and colleagues 27 . This method to analyze bisulfite converted DNA is based on base specific cleavage of nucleotides and subsequent analysis of fragments by MALDI-TOF MS. Methylation values are given for individual CpG sites or multiple CpG sites are analysed together as CpG units depending on the cleavage pattern. The reproducibility of the methodology is remarkably consistent and the methylation values highly accurate 27,28 . According to the results, we found that the average methylation values of circular mtDNA without treatment of BamHI were from 3.43% to 5.93% for CpG106~CpG109, 2.79% to 27.79% for CpG225~CpG243, 1.21% to 12.57% for CpG360~CpG371, while that of linear mtDNA treated with BamHI before bisulfite conversion were from 1.07% to 3.07% for CpG106~CpG109, 1.86% to 6.86% for CpG225~CpG243, 0.57% to 2.29% for CpG360~CpG371, which were significantly reduced ( Table 4). These results further confirm the pyrosequencing results, indicating that the methylation values might be overestimated under the condition of circular mtDNA without cleavage before bisulfite conversion. For the same CpG sites, the methylation values obtained by EpiTYPER were little higher than that of pyrosequencing, e,g, 1.07% to 3.07% vs 0.25% to 0.33% for CpG106~CpG109. It was reported that fully methylated or completely unmethylated samples show consistently a small decrease or increase in methylation, respectively 28  There was no significant difference between two groups for any of the four sequences (Fig. 1). And the pyrosequencing success rate of circular mtDNA group in saliva samples was significantly higher than that in blood samples. Even for the circular mtDNA in saliva samples that was unable to be pyroseqenced, none were caused by incomplete bisulfite conversion ( Table 2). The above results showed that the circular structure of mtDNA has little effect on bisulfite conversion for saliva samples.  For saliva samples, the average methylation values of four sequences were 1.17% ± 0.25%, 1.79% ± 0.35%, 0.47% ± 0.17%, 1.06% ± 0.28% in circular mtDNA group, and 1.38% ± 0.21%, 2.14% ± 0.29%, 0.69% ± 0.18%, 0.93% ± 0.22% in linear mtDNA group, respectively. Wilcoxon signed-rank test results showed that no significant differences exist between the two groups (Fig. 2). For the same samples measured by both methods and for methylation data was obtained, paired comparison was conducted and no statistical difference was found in most of the CpG sites except for CpG8, CpG39, CpG109 sites (Table 3). These results further confirmed that the circular structure of mtDNA in saliva samples has little effect on bisulfite conversion and the subsequent methylation quantification by pyrosequencing. From these results, we understood our previous study that the methylation level of 6 CpG sites in D-loop regions is 2~34% in blood samples but nearly 0% in saliva samples. Because we didn't perform opening-loop treatment before bisulfite conversion, the blood samples underwent an incomplete conversion while the saliva samples were not affected.
The inconsistency between blood samples and saliva samples prompt us to study other tissues. We collected 7 lung tissues from autopsies of forensic cases. The death causes of these 7 individuals were not related with lung diseases. According to the results, among the 6 CpG sites we detected, the methylation values of CpG4 and CpG5 is significantly higher in mtDNA without BamHI treatment than with BamHI treatment. Although no statistical significance was found in other 4 CpG sites, the CpG mean methylation values of the 6 CpG sites showed obvious difference between loop mtDNA group and linear mtDNA group ( Table 5). The methylation value in loop mtDNA group was higher than that of linear group. These results further verified that the loop structure of mtDNA might affect efficient bisulfite conversion.
However, why saliva samples showed different results? Firstly, the saliva DNA amount we used to perform bisulfite conversion was 500 ng, while that of blood and tissue samples was 1000 ng. The excess concentration of DNA in the reaction may interfere with the ratio of bisulfite to DNA. Secondly, we compared the DNA extraction methods of blood, tissues, and saliva samples. We used QIAamp DNA kit to extract DNA from blood and tissue samples, while for saliva samples, silica and guanidinium isothiocyanate method was performed. We noticed that the pH of TE buffer used to purify blood and tissue DNA was 9.0, while that used for saliva cell DNA was 8.0. In the first step of bisulfite conversion, sulfonation, bisulfite is added to cytosine and the reaction is controlled by pH. Low pH induces formation of cytosine sulfonate while high pH reverses the reaction, with the recommended pH being around 5 to keep the process favouring efficient conversion of cytosine to uracil 29 . Therefore, relatively lower pH of saliva DNA might be one of the reasons leading to more efficient bisulfite conversion than blood and tissue samples. For efficient bisulfite conversion, the initial denaturation step is the most critical step as deamination of cytosine to uracil requires DNA to be single-stranded 29 . Incomplete denaturation or reannealing of DNA can result in double-stranded DNA (dsDNA) preventing efficient conversion of unmethylated cytosines. Increased dsDNA formation can also be due to an excess concentration of DNA in the reaction which interferes with the ratio of bisulfite to DNA and also increases the pH 29 . Therefore, if the loop structure of mtDNA is not opened completely, it will be difficult to achieve complete denaturation and form single strand. As a result, the bisulfite conversion efficiency will be reduced. All in together, no matter what sample types or which kind of DNA extractions we used for detect mtDNA methylation, it is best to fragment DNA by using restriction endonucleases that will not cut the sequence of interest to insure the bisulfite conversion efficiency.

Human mitochondrial genome DNA methylation profiling. Applying the developed methods of
BamHI treatment and bisulfite pyrosequencing, we detected mitochondrial genome methylation of blood samples and saliva samples from 14 unrelated individuals. This detection system covers most part of mtDNA including D-loop region, 12 S rRNA, 16 S rRNA, ND1, COXI, ND3, ND4, ND5, CYTB, which can measure 83 CpG sites including 16 CpG sites in the D-loop region, 21 CpG sites in two rRNA gene regions, and 46 CpG sites in six protein coding regions. The results showed that, regardless of sample types, other than two CpG sites within D-loop region, the methylation levels of the rest of the 81 CpG sites were all less than 3%. The mean methylation levels at CpG16 site (nt544) and CpG429 site (nt16411) within D-loop region were relatively high and obviously prominent, which were both about 5% in blood samples and saliva samples. (Supplementary Table 3). The mean methylation values of nine regions were all below 2% (Fig. 3). Other than 12 S rRNA (P = 0.007) and ND1 (P = 0.001), there were no significant differences between different tissues at most detected regions (Fig. 3).
Among the nine detected regions, some of them have never been reported, and it is the first time that methylation status has been obtained for ND1, COX1, ND3, ND4, ND5, and CYTB. However, we cannot rule out the possibility that other CpG sites we didn't analyze may be methylated at high level. But at least, COII and ATP6 have been reported to be unmethylated by Hong et al. 20 . Because bisulfite pyrosequencing cannot measure the methylation of non-CpG sites, we also cannot exclude the possibility that there exists methylation outside CpG in mtDNA. Bellizzi et al. 15 has reported that mtDNA methylation of mouse and human blood or cultured cells particularly occurred within non-CpG nucleotides. However, Hong et al. 20 analyzed published genome-wide bisulfite sequencing data sets and the results indicated the lack of cytosine methylation in mtDNA either at CpG site or non-CpG sites, which was also be confirmed by next-generation sequencing. Several studies documented cytosine methylation in mtDNA by using affinity-based methods such as immunoprecipitation 14 or ELISA 12,13 , did not provide any information on the sites or frequency of CpG methylation of the targets. Our data provided the accurate methylation values at each detected CpG site across most regions of human mtDNA. Except two independent CpG sites in D-loop region, the methylation degree of the remaining 81 CpG sites are all less than 3%. The mean methylation levels at CpG16 site (nt544) and CpG429 site (nt16411) within D-loop region were relatively high and obviously prominent, which were both about 5% in blood samples and saliva samples. It has been reported that the interactions of mitochondria DNA-methyltransferase1 (mtDNMT1) with mtDNA appear to be CpG dependent and particularly evident in the D-loop region 14 . D-loop region carries the mitochondrial H strand origin of replication (O H ). Therefore, higher methylation status at this region might be related with the mtDNA replication, which is need further studies to confirm.
In addition, the sensitivity of pyrosequencing technology for analyzing CpG methylation is 5% and the default value of passed quality of bisulfite treatment is 4.5% at non-CpG cytosine. If the methylation value is close to the threshold, we cannot distinguish the methylation signal from incomplete bisulfite conversion signal. Therefore, the methylation value less than 5% is nonsense and can be considered as non-methylated, and only if the methylation value is lager than 5% under passed quality pyrograph, the methylation result is reliable. With respect to our results, for linear mtDNA, the average CpG methylation degree of 9 detected regions were all less than 5%,  Table 5. Comparison of methylation levels of the same lung DNA sample with or without BamHI treatment detected by bisulfite pyrosequencing.
indicating that there was no CpG methylation at these regions of mtDNA. Furthermore, understanding of the biological impact of DNA methylation requires the critical information of the DNA methylation pattern, including the methylation status of contiguous sites 20 . With regard to the whole regions, the mean methylation values of nine regions are all less than 2%. Therefore it seems that such low level of mtDNA methylation would have limited or absent functional significance in the control of mitochondrial gene expression, as previously suggested by Shmookler 9 . Why is DNA methylation in mtDNA very low or even absent? In vertebrate nDNA, a methyl group is added to the 5′ position of the base cytosine to generate 5 mC by the action of DNMT. In this process, S-adenosylmethionine (SAM) is required as a methyl donor. Recently some studies reported that mitochondria contain the machinery required to epigenetically modify mtDNA 11,14,17,[30][31][32][33] . However there is some discordancy with regard to these results. Shock et al. 14 identified the mitochondrial targeting sequence for DNMT1 in mouse and human, which drives the translocation of mtDNMT1 to the mitochondria, but did not observe the presence of DNMT3a or DNMT3b in mouse and human cells. However, a study by Wong et al. 17 demonstrated the presence of DNMT3a in mitochondria of adult mouse CNS, skeletal muscle, and testes and human cerebral cortex, while DNMT1 was not detected in adult mouse CNS or skeletal muscle mitochondria. Although it has been reported that cytosolic SAM can be transported into mitochondria by the SAM carrier (SAMC), the exact role of SAM into mitochondria remains poorly understood 11,33 . Therefore, whether there is methylation machinery in mitochondria still needs to be reevaluated.
In conclusion, we are the first to show that the circular structure of mtDNA has a great impact on bisulfite conversion efficiency, and thereby affects subsequent pyrosequencing success rate. Even though the sequencing is passed, the quantitative methylation value is significantly higher than the actual value, leading to overestimation of mtDNA methylation. This effect caused by loop structure occurs in blood samples but not in saliva samples. Using BamHI treatment and bisulfite pyrosequencing, human mtDNA methylation profiling of 83 CpG sites covering nine regions (D-loop, 12 S rRNA, 16 S rRNA, ND1, COI, ND3, ND4, ND5, CYTB) was mapped and the average methylation levels of all nine regions are below 2%, indicating CpG methylation in human mtDNA is a rare event at most regions.

Materials and Methods
Samples. According to the principle of informed consent, peripheral blood samples were collected from 154 healthy unrelated individuals, and saliva samples were collected from 24 healthy unrelated individuals. For collection of saliva, the volunteers were asked to rinse mouth with mouthwash and 3 mL saliva was collected to a sterile tube. All the participants signed informed consent. In addition, 7 lung tissues were collected from the autopsies of forensic cases under the consent of their most direct relatives. The death causes of these 7 individuals were all not related with lung diseases. This study passed the ethical review and approved by Hebei Medical University Biomedical Ethics Committee. All of the methods were carried out in accordance with the approved guidelines.
Bioinformatic analysis and primer design of mitochondrial genome CpG sites. Bioinformatic analysis of mitochondrial genome infers a total of 435 CpG sites unevenly distributed across the whole mtDNA, with the highest CpG density in the light chain replication origin (O L ) (5721~5798 bp), and the minimum density in tRNA genes (Supplementary Table 1). Twenty pairs of primers for amplifying specific mtDNA segments and the corresponding sequencing primers were designed using Assay Design Software. These 20 specific products named MT1~MT20 cover both D-loop region and coding region involving 83 CpG sites. The basic information of 20 sequences is shown in Table 1, and the primer sequences are shown in Supplementary Table 2. DNA extraction and quantification. Whole genomic DNA (including mtDNA) was extracted from blood samples and lung tissues using QIAamp DNA Mini and Blood Mini kit (Qiagen, Germany) according to the manual book. DNA extraction form saliva samples was performed by silica and guanidinium isothiocyanate method (Taizhou Stellar Biotechnology Services Ltd, China). DNA concentration and purity were analyzed by ND1000 nucleic acid protein quantitative analyzer. The average concentration of DNA in blood samples and saliva samples were 171 ng/μl and 94 ng/μl, and the average OD260/280 for both samples were 1.90, 1.76 respectively. DNA integrity was verified by 1% agarose gel electrophoresis using randomly selected 60 blood samples and all of the saliva samples.
BamHI treatment. BamHI specifically recognizes GGATCC and cuts at the position of G/G. There is only one BamHI restriction site in mtDNA at nt14258~nt14263. So after BamHI treatment, the loop of mtDNA can be opened and form a complete linear DNA. DNA was incubated with BamHI at 37 °C for 4 hours in a 20 μl reaction system comprising DNA 1 μg, 10× Buffer 2 μl, BamHI 2 μl, to get complete cutting.
Bisulfite conversion, PCR and pyrosequencing. Bisulfite conversion was performed using the EZ DNA Methylation-Gold Kit (Zymo Research, USA) according to the manufacturer's protocol. The converted DNA was amplified to get specific mtDNA fragments using specific primers (Supplementary Table 2) and GoTaq Green Master Mix (Promega, USA). The total reaction volume was 50 μl containing GoTaq Green Master Mix 25 μl, nuclease free water 22 μl, 10 μM sense primer and anti-sense primer 1 μl, and converted genomic DNA 30 ng. The cycle conditions were an initial denaturation at 95 °C for 5 min, 48 cycles of 95 °C for 15 s, 51~57 °C for 30 s, 72 °C for 15 s, and a final extension at 72 °C for 5 min. PCR products were subjected to gel electrophoresis on 2% agarose. If we obtained a specific fragment, then the PCR products could be analyzed by pyrosequencing as follows.
The single strand of PCR products coupled with biotin was purified by mixing 37 μl PCR products with 40 μl binding buffer and 3 μl avidin beads for 30 min with gentle shaker at room temperature, followed by washing with 70% ethanol for 5 s, denaturing with denaturation buffer for 15 s, and washing with washing buffer for 15 s. The purified single strand was annealed with 1.6 μl sequencing primer and 38.4 μl annealing buffer at 80 °C for 5 min. The sequences of the sequencing primers were shown in Supplementary Table 2. After setting up the inner control for testing the efficiency of bisulfite conversion, the pyrosequncing was conducted using PyroMark ID and the quantitative methylation level of each CpG site were obtained. All of the above methods were carried out in accordance with the approved guidelines.

SEQUENOM EpiTYPRE MassARRAY.
To validate the results of pyrosequencing, the mtDNA methylation status of fourteen blood samples was reanalyzed using SEQUENOM MassARRAY platform (CapitalBio, Beijing, China). PCR primers were designed with Methprimer (http://epidesigner.com). For each reverse primer, an additional T7 promoter tag for in vivo transcription was added, as well as a 10-mer tag on the forward primer to adjust for melting temperature differences. Three pairs of primers were designed and used to amplify 3 fragments of mtDNA, from which 31 CpG sites were detected totally. The detailed information was listed in Supplementary  Table 4. The bisulfite converted DNA was amplified using PCR Accesory set (Sequenom, San Diego, CA) and the amplification products were cleaned with shrimp alkaline phosphatase using Mass CLEAVE kit (Sequenom, San Diego, CA) to dephosphorylate all unincorporated nucleotides. After purification in vitro transcription of the reverse strand and base-specific (C or T) cleavage by RNase A was performed. The cleavage products, which differ in mass and length were spotted onto a MALDI matrix-containing SpectroCHIP and subjected to MALDI-TOF MS. The mass spectra was collected by MassARRAY Spectrometer and analyzed by EpiTYPER software version 1.0 (Sequenom, San Diego, CA). Statistical analysis. Statistical analyses were performed using the SPSS13.0 statistical software. Paired t test, variance analysis, Wilcoxon test were adopted, and a two-sided P < 0.05 was considered statistically significant.