## 2 Materials and Methods
### Bioinformatic Analysis
The combined approach enabled identification of haplotype-specific methylation signatures that would be undetectable by analyzing genomic variants or methylation patterns alone, providing insights into potential functional consequences of mitochondrial sequence variations. Visualization of these relationships was accomplished through heatmaps and profile plots comparing methylation patterns between haplotype clusters, providing insights into potential functional relationships between sequence variants and epigenetic modifications in the mitochondrial genome. Additionally, several analytical data analysis strategies were implemented to determine if the mtDNA methylation phenomena observed was an artifact or a result from genuine 5 meC methylation. First, to investigate the relationship between mitochondrial DNA (mtDNA) methylation and sequencing coverage obtained from adaptive sampling Nanopore long-read sequencing, we developed a custom Python analysis pipeline. The pipeline utilized aligned reads stored in BAM format as input. Initially, the mitochondrial chromosome was identified within the BAM header, and the presence of standard methylation call tags (e.g., MM, ML, YM) was automatically detected. Reads mapping to the mtDNA were processed using the pysam library (version dependent on environment). Reads were filtered based on a minimum mapping quality score (MAPQ \(\geq\) 20) and exclusion of secondary or supplementary alignments. For each filtered read, per-base coverage information was recorded. Methylation status (methylated or unmethylated) for 5 meC sites (or other relevant contexts depending on tag parsing logic) was extracted by parsing the appropriate BAM tags, converting relative read coordinates to genomic positions. Data were aggregated across all reads to calculate the total coverage and the frequency of methylation (number of methylated reads/total reads assessed) at each genomic position. Only positions meeting a minimum coverage threshold (e.g., \(\geq\)10 reads) were retained for downstream analysis. The relationship between per-position coverage and methylation frequency was assessed using Pearson correlation, Spearman rank correlation, and linear regression analysis (scipy.stats, sklearn.linear_model). To evaluate relationships at different genomic scales, position-level data were further aggregated into non-overlapping genomic bins (e.g., 50, 100, 250, 500 bp), calculating the average coverage and average methylation frequency per bin. Correlation analyses were repeated on these binned data. Potential coverage biases were assessed by stratifying positions into quantiles based on coverage and comparing methylation frequencies, and conversely, stratifying by methylation frequency categories and comparing average coverage using bar plots, box plots, and one-way ANOVA (scipy.stats). _Methylation Pattern Analysis_: Methylation status was assessed by extracting base modification information from BAM tags (e.g., MM, ML, XM, YM), typically generated by basecallers or post-processing tools like Nanopolish, Megalodon, or Guppy+modkit. For reads aligned to the mtDNA reference, methylation frequencies were calculated within each 100 bp bin.