## 2 Materials and Methods
### Bioinformatic Analysis  
Oxford Nanopore bioinformatic workflow Nextflow was used for super-accuracy model re-basecalling (wf-bascell, Dorado v0.7.0, SUP model) and human genome alignment (wf-alignment, Ensembl reference Hs.GRCh38 (release 110)) of the two produced sets of POD5 files for the libraries. Python code, using pysam and mathplotlib, and NumPy was used to generate circus-like visualization displaying mtDNA coverage depth throughout the mitochondrial genome, with integrated gene annotations and statistical highlighting of over-enriched regions relative to the nuclear genome. Statistical significance for enrichment was determined by calculating a threshold based on 2\(\times\) the median enrichment across all mtDNA genes, where enrichment was defined as the ratio of a gene's mean coverage to the average nuclear genome coverage sampled across multiple chromosomes. Mitochondrial DNA methylation analysis was performed through a multi-step pipeline for analysis of the two WGS and adaptive sampling BAM files (obtained above). Cytosine modifications were then extracted with a custom Python 3.10 notebook that uses pysam v0.22.0, NumPy v1.26, SciPy v1.11 ('scipy.stats.spearman'), Matplotlib v3.8 and tqdm v4.66. The algorithm (i) retains only C\(\rightarrow\)5-mC or C\(\rightarrow\)5-hmC calls whose 'ML' posterior probability \(\geq\) 0.80, (ii) maps read coordinates to reference positions, (iii) classifies a site as non-CpG when the reference dinucleotide is C H (H = A, C, T), and (iv) bins methylated versus total CH calls in 100 bp windows to calculate % 5-modC. Mean read-depth and window-level methylation are plotted along the 16.6 kb mitogenome, and a methylation-versus-coverage scatter-plot--annotated with Spearman \(\rho\)--confirms that the residual signal is independent of sequencing depth. Heteroplasmic variant analysis was conducted by examining the mitochondrial genome for positions showing significant alternate allele frequencies (\(\geq\)5%). The analysis pipeline processed BAM files using pysam, applying quality filters (mapping quality \(\geq\) 20, base quality \(\geq\) 20) to identify reliable variants. The pipeline built a read-variant matrix to analyze the distribution of heteroplasmic variants across individual reads. We employed K-means clustering with silhouette score optimization to identify distinct mitochondrial subpopulations (K-means k = 3 chosen by silhouette score = 0.71; first two PCs explaining 87% variance). Principal component analysis (PCA) was used to visualize the clustering of reads based on their variant profiles. For each identified cluster, we calculated allele frequencies at heteroplasmic sites and created consensus sequences. Reads covering heteroplasmic sites were then extracted to reconstruct potential haplotypes, with reads grouped by their shared nucleotide patterns at variant positions. In samples with insufficient heteroplasmic sites for standard haplotype reconstruction, methylation patterns were employed as alternative clustering features using K-means and hierarchical clustering algorithms. Integration of epigenetic modifications with haplotype data was performed to identify potential correlations between sequence variations and methylation patterns. For each identified haplotype, methylation frequencies were calculated at each position and compared across haplotypes to identify differentially methylated regions (DMRs). DMRs statistical significance was assessed through custom algorithms designed specifically for linking epigenetic modifications to sequence variations, by comparing methylation differences between haplotypes, with positions showing \(\geq\)30% methylation difference highlighted. These algorithms accounted for the circular nature of mtDNA and potential strand-specific methylation patterns.