Single-cell multi-omic analysis of mitochondrial mutational mosaicism and dynamics

Domain: biology
Published: 2026-03-16
Editorial outcome: minor_revision

Abstract
Mitochondrial DNA (mtDNA) mutations occur more frequently than nuclear mutations and are associated with various diseases. While single-cell sequencing enables mtDNA variant heteroplasmy analysis, a holistic view of mtDNA mutational landscapes in individual cells has remained limited. Here, we leverage mitochondrial single-cell ATAC-seq and mtDNA-hypermutated POLGD274A knock-in HEK293 cell lines to introduce two metrics—single-cell mtDNA mutations per million base pairs (scmtMPM) and heteroplasmy-weighted mitochondrial local constraint scores (scwMSS)—to capture cellular mutational loads and somatic mosaicism. We demonstrate that individual POLGD274A cells exhibit complex mutational landscapes, with pathogenic mutations and truncating variants only present at subthreshold levels, indicative of their negative selection. In human healthy donors and mitochondriopathy patients, we identify constrained mutations in complex I, highlighting previously unrecognized mtDNA mutational landscape heterogeneity present on the single-cell level. Overall, scmtMPM and scwMSS provide a framework to investigate fundamental properties of mitochondrial genetics, disease, and somatic mosaicism.

Introduction
Mitochondria are essential for cellular metabolism and contain a multi-copy genome with hundreds to thousands of copies per cell1,2. In humans, the 16.6 kb circular double-stranded mitochondrial DNA (mtDNA) encodes 13 protein subunits, contributing to four complexes of the respiratory chain, alongside its own translational machinery, comprising two ribosomal RNAs (rRNAs) and 22 transfer RNAs (tRNAs). Notably, somatic mtDNA mutations occur more frequently than in their nuclear counterpart, with an up to >200× fold elevated rate in human blood3. Bulk-sequencing has demonstrated an age-related increase in somatic mtDNA variants4,5, established multiple associations with cancers6–8, and other human phenotypes1. However, establishing genotype-phenotype relationships at the cellular level remains challenging due to heteroplasmy (the presence of wild-type and mutant mtDNA) with variable variant allele frequencies (VAFs), and cell-type- and state-dependent mtDNA copy number. Additionally, relaxed replication, vegetative segregation, and natural turnover of mitochondria collectively make cellular mtDNA VAFs more dynamic than for nuclear DNA variants2.

The pathogenicity of mtDNA variants is best established in the context of maternally inherited mitochondriopathies9. While typically multi-systemic, cell-type-specific effects may extend to red blood cell production, and selection against pathogenic mtDNA has been described in T cells10–13. Defects in the murine mtDNA polymerase lead to a ‘mutator’ phenotype with increased levels of mtDNA variants that result in the premature onset of aging-related features14,15. However, fundamental challenges remain in our understanding of the functional impact of mtDNA variants. While individual pathogenic single-nucleotide variants (SNVs) or large-scale deletions are readily definable, the additivity of constantly arising somatic variants is harder to account for. Depending on the specific gene affected and its VAFs, these variants may exert variable cell-type-specific effects. Traditionally, cells are thought to tolerate 60–80% of mutated mtDNA before reaching the biochemical threshold at which metabolic defects manifest themselves1,16; however, recent genetic evidence suggests that phenotypes may already manifest at lower VAFs11,17.

Recent advancements in single-cell sequencing have enabled massively parallel whole mitochondrial genome sequencing18, and have demonstrated an unprecedented level of variation in mtDNA variants19 and copy number20 across thousands of cells. For example, mitochondrial single-cell ATAC-seq (mtscATAC-seq) enables the concomitant sequencing of mitochondrial genomes and chromatin accessibility19,21. Moreover, a recently introduced mitochondrial genome constraint model provides a means of interpreting genetic variation underlying human phenotypes22 across genomic databases such as the Genome Aggregation Database (gnomAD). Here, we leverage these advances to characterize cell lines with proofreading-deficient mtDNA polymerases and to benchmark the sensitivity of these approaches for detecting and quantifying mitochondrial genetic variation. We demonstrate an even more substantial increase in mutational load than previously appreciated via bulk-sequencing and introduce a computational framework to quantify and interpret the genome-wide mtDNA mutational burden. We apply this framework to peripheral blood mononuclear cells (PBMCs) from human donors and mitochondriopathy patients to aid the identification of pathogenic mtDNA variants. Our widely applicable framework opens avenues for single-cell-resolved study of mitochondrial genetics and somatic mosaicism in primary human cells and tissues.

Results
Elevated single-cell mtDNA burden in POLGD274A cells
To assess the sensitivity of mtscATAC-seq to quantify increases in relative mitochondrial genome-wide mutational burden at the single-cell level, we harnessed exonuclease-deficient DNA polymerase γ (POLGD274A) knock-in HEK293 cell lines with an elevated mutational rate and hypermutated mitochondrial genomes23–25. Both KI36 and KIA2 lines showed metabolic impairment compared to CTRL (Supplementary Fig. 1a). The control and POLGD274A cell lines were jointly profiled using mtscATAC-seq with cell hashing to distinguish individual lines and detect mtDNA variants alongside accessible chromatin profiles (Fig. 1a; Supplementary Fig. 1b, c; “Methods”). Following deep sequencing, quality control metrics were comparable among all three cell lines, which clustered distinctly, indicative of heterogeneous chromatin accessibility (Fig. 1b; Supplementary Fig. 1d). Of note, single-cell mtDNA depth was 3–4-fold higher in both POLGD274A knock-in lines (Fig. 1c; Supplementary Fig. 1d; median raw mtDNA depth CTRL: 45, KI36: 130, and KIA2: 214). Given that KIA2 cells also exhibited higher nuclear read counts (median single-cell ATAC-seq read counts CTRL: 18,398, KI36: 17,113, and KIA2: 41,700) that generally correlated with mtDNA depth across all three lines (Spearman correlation > 0.5), the mtDNA depth of individual cells was normalized to their respective nuclear ATAC read count (Fig. 1d; Supplementary Fig. 1d, e). However, the normalized single-cell mtDNA depth in both POLGD274A knock-in lines remained significantly higher than in CTRL (Fig. 1e, f), suggesting higher absolute mtDNA content, likely reflecting a compensatory response in copy number to an elevated mutational load25,26. Notably, in POLGD274A cells, the pseudobulk rolling-mean mtDNA coverage appeared lower in the minor arc, the smaller segment of the genome between the replication origins of the heavy strand (oriH) and the light strand (oriL). This uneven distribution is consistent with the findings that POLG 3′−5′ exonuclease deficiency lead to the accumulation of linear mtDNA fragments in the major arc (Fig. 1g)23,24.

Using the mgatk19 pipeline, 620 mtDNA variants were identified across CTRL cells. In the KI36 and KIA2 cells, a ~15-fold increase in mtDNA variants were identified, 9656 and 11407, respectively (Fig. 1h; Supplementary Data 1). While this higher number of mtDNA variants following disruption of POLG proofreading deficiency is expected, the overall magnitude is more pronounced than previously described, where only hundreds of somatic variants were identified by bulk-sequencing14,27,28. Both control and POLGD274A cells displayed a C>T transition-dominated mutational signature, indicating replication error as the primary driver (Fig. 1i)5,27,29. Notably, relaxing the default thresholds in the mgatk variant calling resulted in a substantial increase in detected variants (Fig. 1h, blue dotted box; “Methods”), which down to a particular strand concordance still showed the expected mutational signature (Supplementary Fig. 1f). Presuming similar background levels compared to control, these findings suggest an even higher number of bona fide somatic mtDNA variants in POLGD274A cells. Together, these observations demonstrate that mtscATAC-seq is a sensitive tool for detecting elevations in mtDNA mutational rates and the resulting increase in mtDNA variants.

Adaptive gene regulatory responses in POLGD274A cells
To dissect nuclear–mitochondrial interactions in response to elevated mtDNA mutational burden, we performed differentially accessible gene (DAG) analysis of POLGD274A vs. CTRL cells. For each POLGD274A line, we identified 827 and 2319 upregulated DAGs, and 1198 and 3113 downregulated DAGs for KI36 and KIA2, respectively (Supplementary Fig. 2a; Supplementary Data 2, 3). Among these, 243 upregulated and 277 downregulated DAGs were shared between the KI36 and KIA2 lines. Upregulated genes included ones in extracellular matrix remodeling and mechanotransduction30 (e.g., CYR61, TGFBI, ADAMTS1), protocadherins (PCDH gene family), and multiple zinc-finger (ZNF) transcription factors31, indicating altered cell–matrix interactions and epigenetic rewiring in POLGD274A cells. Interestingly, mitochondrial genome maintenance exonuclease 1 (MGME1), a mitochondrial exonuclease essential for mtDNA maintenance and transcription associated with mitochondrial disease32, was among the downregulated genes in POLGD274A cells (Supplementary Fig. 2b). MGME1 deficiency23,33 has been shown to lead to accumulating linear mtDNA molecules, and Mgme1-knockout in mice elicited a strong inflammatory response via the cGAS–STING pathway34. Additionally, several genes involved in cellular metabolism (e.g., SLC25A31, SLC25A53, OTC) and respiratory chain components/assembly35 (e.g., COX8A, TIMMDC1, SAMM50) were downregulated in POLGD274A lines, suggesting a global attenuation of mitochondrial biogenesis and OXPHOS function. Together, these findings demonstrate that mtscATAC-seq captures nuclear-mitochondrial interactions linked to mitochondrial dysfunction.

VAF-dependent mtDNA variants dynamics in POLGD274A cells
Given the multi-copy nature of the mitochondrial genome, mtDNA variants may be present only in a fraction of the copies within a cell or population, a state known as heteroplasmy. As mtscATAC-seq yields a high mtDNA depth per cell (up to 400x and 800x in control and POLGD274A knock-in cells, respectively), we determined the number of variants with a VAF exceeding a specified threshold. First, we detected a median of 35 high-confident mtDNA variants per CTRL cell and 904 and 1487 in the two POLGD274A cell lines regardless of their VAFs (Fig. 2a). The total number of mtDNA variants per cell primarily correlated with mtDNA sequencing depth, particularly at lower VAF levels (Supplementary Fig. 3). However, the mtDNA-depth-dependent increase in mtDNA variants detected per cell begun to plateau for variants >3% VAF and higher and became mostly independent of depth >10% VAF (Supplementary Fig. 3). This aligns with the expectation that the total detected number of low VAF variants will be a function of recovered mtDNA reads. In contrast, for higher VAF variants, their confident detection does not further benefit from deeper coverage. For further analysis, we opted to normalize the total number of mtDNA variants per cell by cellular mtDNA sequencing depth (Fig. 2b; “Methods”), after which the two POLGD274A lines still showed an increase in the depth-normalized number of mtDNA variants per cell (median: CTRL: 0.769, KI36: 7.01, and KIA2: 6.97).

To refine the mutational analysis, we only considered ‘shared variants’ between the control and POLGD274A lines, revealing 593 mtDNA parental variants. Notably, KIA2 showed a higher number of shared variants with the parental CTRL line (n=586), compared to KI36 (n=407), likely indicative of their more direct clonal relationship (Fig. 2c). In contrast, variants unique to KIA2 and KI36 cells were classified as “POLGD274A-line specific” variants, which likely arose ‘de novo’, though it is possible that these were already present at very low frequency in the parental line. Filtering cell line-specific variants by their VAFs notably showed a distinct distribution, with the number of low-frequency mtDNA variants, in particular, showing a 50-to 100-fold increase in POLGD274A compared to control cells (Fig. 2d, e). This is consistent with constantly newly arising variants initially only being present in a single mtDNA copy, thus at low VAFs. Over time, relaxed replication and genetic drift can increase their prevalence across additional mtDNA copies. As such, the total number of detected mtDNA variants per cell progressively decreases when higher VAF thresholds are applied (Fig. 2d, e). Along these lines, parental variants displayed higher pseudobulk VAF levels than “line-specific” variants (Fig. 2f, g). These observations align with anticipated heteroplasmy dynamics of somatic events, with higher VAF levels requiring further genetic drift and vegetative segregation36. Together, these analyses demonstrate that mtscATAC-seq sensitively detects de novo somatic mtDNA variants in single cells, even at low allelic frequencies.

To further examine the VAF distribution of the most variable variants, we categorized them into synonymous, missense, truncating, and pathogenic mtDNA variants previously reported in MITOMAP37 (Fig. 3a–d; “Methods”). Synonymous and missense mtDNA variants showed broad VAF distributions, with a substantial fraction of variants reaching higher VAF (25–50%) in single cells. Collectively, these variants exhibited a wide range of Shannon entropy values, reflecting diverse heteroplasmy VAF patterns (Fig. 3b, c). In contrast, both pathogenic and truncating mtDNA variants exhibited highly skewed VAF distributions, with nearly all remaining at zero or low allelic frequencies and well below thresholds typically associated with mitochondrial dysfunction of pathogenic variants (e.g., 60–80%)16 (Fig. 3b, c). Moreover, out of the 130 reported pathogenic variants in MITOMAP (Supplementary Data 4), 43 and 41 variants were detected in KI36 and KIA2, respectively. These deleterious variants also exhibited consistently lower entropy values, suggesting restricted heteroplasmy dynamics at sub-threshold levels, indicative of their negative selection during continuous culture. Pairwise assessments of VAF distributions of truncating and pathogenic variants further revealed that many pathogenic variant pairs showed a dense cluster of cells with zero VAFs for both variants and appear only infrequently concomitantly detected in the same cells (Supplementary Fig. 4a). In contrast, truncating variants were more frequently observed together in individual cells, including at higher but mostly <60% VAF, resulting in a greater number of double-positive cells (Supplementary Fig. 4b). In summary, POLGD274A cells show complex mtDNA variant patterns, with truncating and pathogenic mtDNA variants being relatively depleted at higher VAF levels. These variants remain below classical biochemical thresholds in individual cells, suggesting that a complex and widespread mutational burden, as opposed to a single or few deleterious variants, may underlie POLG-mutation associated mitochondrial defects (Supplementary Fig. 1a). While mtscATAC-seq does not directly establish genotype-phenotype associations at the metabolic level, these findings provide important insights into POLG mutator genetics and motivate the systematic investigation of how a cumulative mutational burden, rather than single variants, impact mitochondrial, metabolic, and cellular function.

Single-cell mtDNA mutational burden in POLGD274A cells
To date, most single-cell resolved mtDNA genotyping approaches have primarily focused on individual mtDNA variants rather than considering the complete variant call set per cell, which may further be highly variable across a population of cells11,19. Specifically, individual POLGD274A cells may carry hundreds of mtDNA variants at variable VAFs, and variable pathogenicity distributed throughout the mitochondrial genome, affecting different genes and respiratory complexes (Figs. 3, 4a). Motivated by the heterogeneous VAF distributions of mtDNA variants, we aimed to establish a framework to quantify and infer the functional consequences of the elevated total mtDNA mutational burden of individual POLGD274A cells. Specifically, we introduced two metrics for evaluating the overall single-cell mtDNA mutational burden (scMMB). First, the number of mutations per million base pairs (MPM) is widely recognized for quantifying the mutational burden in the nuclear genomes of cancers38–40. This metric has been adapted for assessing the mitochondrial genome mutational burden in bulk exome sequencing datasets6. By normalizing the number of mutations to sequencing depth, MPM further provides a standardized measure that facilitates comparisons across samples and conditions. Here, we adapt the principle of MPM quantification to single-cell mitochondrial genome sequencing data, which we termed single-cell mtDNA mutations per million base pairs (Fig. 4b, scmtMPM). Because mtDNA depth may vary from cell to cell, scmtMPM provides a normalized assessment of the overall mutational frequency of the mitochondrial genome and its heterogeneity in single cells.

In addition to scmtMPM, a recently introduced constraint model built on large-scale genomic data from gnomAD22,41 compares the observed genetic variations with those expected under neutral evolution to assess the local constraints of the mitochondrial genome. Specifically, a mtDNA local constraint (MLC) score enables the quantification of the local tolerance to base substitutions. The MLC score ranges from 0 to 1, with 1 indicating the most constrained positions and, consequently, the most likely to have a deleterious impact when mutated. Furthermore, the sum of MLC scores (MSS) across a sequenced whole genome enables the inference of the overall implications of mtDNA genome-wide variation, which has further been shown to be significantly associated with human phenotypes and the risk of disease progression22,41,42. Because individual cells may exhibit unique VAF distributions of mtDNA variants2 (Fig. 4a), we further account for their VAFs rather than only considering the presence or absence of a variant. For example, a heteroplasmic variant present at 80% VAF is expected to exert a more substantial phenotypic effect compared to when present at 5% VAF. As such, we propose a single-cell Heteroplasmy-Weighted MSS (scwMSS) metric, which weights the contributions of each variant MLC by its VAF to allow for a refined evaluation of the potential functional impact (Fig. 4b).

For simplicity, we first evaluated scmtMPM and scwMSS scores by functional groups of genes, i.e., encoding for complex I, III, IV, and V, rRNAs, and tRNAs, analogous to prior bulk-sequencing based studies, for example, having described distinct mutational burden across complexes during tumor evolution6. As expected, POLGD274A cells showed substantially elevated scmtMPM and scwMSS scores, which were independent of mtDNA sequencing depth (Fig. 4c, d; Supplementary Fig. 5a–d; median total scmtMPM, CTRL: 0, KI36: 2084, and KIA2: 1796; median total scwMSS, CTRL: 0, KI36: 16.8, and KIA2: 14.7). The presence of cluster-specific mtDNA variants thereby suggests independent subclonal trajectories underlying the acquisition of somatic mitochondrial genetic mosaicism (Fig. 4e).

Generally, the control cells exhibited an overall low mtDNA mutational burden across functional groups, with a few control cells showing elevated scores for complex I and IV (Supplementary Fig. 5b, c). In contrast, POLGD274A cells showed a narrower distribution of scmtMPM and scwMSS scores, which were generally higher in the KIA2 derivative (Fig. 4c, d; Supplementary Fig. 5a–c). Overall, scmtMPM and scwMSS scores were well correlated, though this correlation was more apparent in POLGD274A cells (Supplementary Fig. 5e), suggesting that both measures provide a means of quantifying the cellular genome-wide mtDNA mutational burden.

Further, we assessed correlations in the scmtMPM scores across functional groups and genes across all individual cells (Supplementary Fig. 6a, b). In the control line, correlations were overall weak (r < 0.1), and only a few cells displayed a high burden in more than one respiratory chain complex, suggesting the sporadic independent acquisition of mutations across functional groups and genes. However, in the hypermutated cells, the scmtMPM scores exhibited an elevated scmtMPM score and increased positive correlation (median r = 0.33, interquartile range 0.25–0.55) across functional groups and genes, consistent with a continuous and genome-wide somatic mutational process uniformly affecting the entire mitochondrial genome (Supplementary Fig. 6a, b). Collectively, these data highlight the utility of scmtMPM and scwMSS scores to quantify the genome-wide mtDNA mutational burden at the level of individual cells, complexes, and their interrelationships.

Depth-dependent mtDNA variant discovery and saturation
To evaluate how mtDNA sequencing depth impacts variant detection, we downsampled mtDNA reads to defined depths and quantified the number of detected variants in pseudobulk and in single cells. As expected, the number of total detected variants increased with depth. At shallow depths (<40M reads; Supplementary Fig. 7a), the relationship was approximately linear, with a linear regression model yielding slopes of 17.9, 134.9, and 141 variants per million reads for CTRL, KI36, and KIA2, respectively. This suggests an 8-fold higher discovery rate in KI36 and KIA2 relative to CTRL. At higher sequencing depths (>50M reads, Supplementary Fig. 7b), the rate of discovery plateaued in POLGD274A lines. Next, we evaluated the impact of sequencing depth on mtDNA variant detection at the level of single-cells and observed a similar depth-dependent trend as in pseudobulk (Supplementary Fig. 7a, c, d). Linear regression models on the median number of single-cell mtDNA variants yielded slopes of 1.25, 12.27, and 16.77 variants per million reads for CTRL, KI36, and KIA2, respectively, indicating >10-fold higher discovery rates in KI36 and KIA2 relative to CTRL, consistent with the pseudobulk analysis.

Despite the expected depth-dependent effects on variant number, their heteroplasmy estimates were remarkably stable. Pseudobulk VAFs and variant composition remained concordant across depths (Supplementary Fig. 7e, f). Single-cell mtDNA heteroplasmy quantifications were highly consistent between 50M (corresponding to ~25,000 mtDNA reads per cell) and 200M (~100,000 mtDNA reads per cell) total mitochondrial reads, with most per-variant changes in VAF (ΔVAF) values within ±5% and a global ΔVAF distribution tightly centered at zero (Supplementary Fig. 7g). Likewise, the scmtMPM and scwMSS metrics remained relatively stable across a broad depth range and showed minimal bias even at 20M total mtDNA reads (~10,000 reads per cell; Supplementary Fig. 7h), in contrast to the strongly depth-dependent raw mutation counts. Together, these analyses indicate that while sequencing depth is an important technical factor impacting variant discovery, mtDNA heteroplasmy and mutational burden remain stable once moderate coverage is achieved, and depth is appropriately accounted for.

Galactose amplifies mtDNA copy number without clonal selection
We next examined how POLGD274A cells respond to metabolic stress imposed by galactose43–47. Specifically, galactose is selectively lethal to cells carrying pathogenic mtDNA variants46 and selects against cells with higher large-fragment mtDNA deletion heteroplasmy47. Notably, CTRL cells showed a transient delay in cell growth but proliferated normally by day 3. In contrast, POLGD274A cells showed substantial cell death, with KIA2 cells not tolerating galactose treatment and only a fraction of KI36 cells surviving, indicating a pronounced selective pressure (Supplementary Fig. 8a).

To test whether the selection dynamics occurred at the level of individual variants or the overall mtDNA mutational burden, we performed mtscATAC-seq of CTRL and KI36 cells cultured in 10 mM glucose or galactose for 1 and 3 days, yielding a total of 6165 and 7664 high-quality CTRL and KI36 cells, respectively (Fig. 5a; Supplementary Fig. 9b–d). In line with our baseline characterization of the POLGD274A lines (Fig. 1b), CTRL and KI36 clustered separately (Fig. 5a; Supplementary Fig. 9a, b). Notably, galactose-treated KI36 cells clustered more distinctly in response to galactose, likely reflecting a compensatory response (Fig. 5b). At the level of mtDNA variants, we did not detect notable shifts in pseudo-bulk heteroplasmy levels in either CTRL or KI36, nor did we observe alterations in the overall distribution of single-cell VAFs under galactose compared to glucose conditions (Fig. 5d, e; Supplementary Fig. 9d–f). Moreover, the mutational burden remained comparable between glucose and galactose treatment (Supplementary Fig. 9f). However, there was a significant increase in mtDNA copy number in both CTRL and KI36, with KI36 showing a more pronounced increase (Fig. 5c; Supplementary Fig. 9c; fold change of normalized median mtDNA depth for Gal Day 3 vs. Glu Day 3: 1.3× in CTRL and 2.3× in KI36), indicating that mtDNA amplification represents the primary mode of compensation in this setting47–49.

To investigate correlating changes in chromatin accessibility, we focused on two large subclones marked by specific somatic mtDNA variants and compared their epigenetic state in glucose vs. galactose treatment (Fig. 5f). Differentially accessible gene (DAG) analysis revealed 112 and 63 DAGs in C1 and C3, respectively (Supplementary Data 5, 6). Among these, 10 DAGs were upregulated in galactose treated cells in both C1 and C3 and involved genes in mitochondrial/ER stress and inflammatory signaling pathways (e.g., CCL20, TRIB3, INHBE, ARHGEF2)50–52, and four non-coding RNAs (i.e., DLGAP1-AS2, MIR6738, MIR6864, and MIR6865), likely reflective of metabolic rewiring and adaptation to OXPHOS-enforcing conditions in galactose (Fig. 5g). Interestingly, GDF15, a well-known biomarker of mitochondrial diseases53–55, was also identified as DAG in C1, showing a gradual increase in chromatin accessibility under continuous galactose treatment (Fig. 5h). As strongly pathogenic mtDNA variants were largely absent in this context, galactose appeared to have a predominantly uniform effect across the population rather than selecting for a favorable mtDNA variant profile, with elevation in mtDNA copy number being the primary mode of adaptation.

Heterogeneous mtDNA somatic mosaicism in human blood cells
Maternally inherited mtSNVs and mtDNA somatic variants acquired during aging shape the mutational landscape and have been implicated in human phenotypes, including aging, metabolic disease, neuromuscular disease, cancers, and congenital mitochondrial diseases1,2,7,9,56. Moreover, the co-occurrence and mutual exclusivity of mtSNVs suggest complex interactions among mtDNA variants19,56. However, the dynamics of mtSNVs and the impact of pathogenic mtDNA variants on the mtDNA mutational landscape remain largely unexplored. To address these questions, we first evaluated the applicability of scMMB scores to primary human cells. We analyzed PBMCs from two healthy individuals (5- and 47-year-old) and four patients (29-, 35-, 60-, and 80-year-old) with mitochondrial encephalomyopathy with lactic acidosis and stroke-like episodes (MELAS) syndrome (Fig. 6a). Across all individuals, scmtMPM and scwMSS scores were well correlated and independent of mtDNA depth (Supplementary Fig. 10). Additionally, an age-related increase in the number of detected mtDNA variants was observed across healthy donors and patients (Supplementary Fig. 11a, b), consistent with prior observations3. However, a more heterogeneous distribution of scmtMPM and scwMSS scores in primary cells, compared with cell lines, was observed in a cell-type- and donor-specific manner (Fig. 6b, c; Supplementary Fig. 5).

Specifically, in the pediatric donor, groups of cells with a relatively higher mutational burden were observed across different immune cell types (Fig. 6c; top left, arrows), suggesting the non-uniform acquisition and distribution of somatic variants across cell types. As expected, the immune cells of patients with pathogenic mtDNA variants showed higher scmtMPM and scwMSS scores and a more substantial bimodal distribution (Fig. 6b; red violin plots). As the pathogenic mt.3243A>G mutation undergoes purifying selection in specific T cell subsets and the hematopoietic system throughout the lifespan of patients (Supplementary Fig. 11c)10,11,13, we also reevaluated scMMB metrics after excluding the pathogenic variant. Notably, the removal of the mt.3243A>G mutation abolished the bimodal distribution in the younger patients (29- and 35-year-old). However, this reduction was not observed in the older patients (60- and 80-year-old), suggesting that while the mt.3243A>G mutation may be gradually lost, new mutations emerge with age contributing to elevated scmtMPM and scwMSS scores (Fig. 6b). Additionally, a significant fraction of cells with no detectable mt.3243A>G heteroplasmy, showed relatively higher scmtMPM scores (Fig. 6d, e; arrows) compared to cells even with low mt.3243A>G VAFs, suggesting the presence of a highly pathogenic congenital mtDNA variant to potentially exert selective pressures to accumulate additional somatic mtDNA variants. Notably, reevaluation after ‘removing’ mt.3243A>G showed that the scmtMPM metric remains low in cells with detectable levels of the variant, regardless of mt.3243A>G VAFs, mtDNA depth, and cell type (Fig. 6e, f; Supplementary Fig. 11d). Together, scmtMPM and scwMSS scores enable the single-cell resolved evaluation of patterns of pathogenic mtDNA variants and demonstrate their impact on somatic mosaicism in primary human cells.

Increased complex-specific mtDNA mutational burden across lineages
To further investigate the nature of elevated scMMB scores, we evaluated whether they are attributable to a single highly constrained variant or to the cumulative effect of multiple mtDNA mutations. For example, a single highly constrained variant would over proportionally contribute to elevated scMMB scores. To assess this, we deliberately removed the variant with the highest VAF heteroplasmy, after which scMMB scores decreased significantly, suggesting that a single somatic variant profoundly affected scMMB scores (Fig. 7a; yellow violin plots). Next, we stratified scores based on gene and functional groups to further resolve the underlying heterogeneity of scmtMPM and scwMSS in cells with pathogenic mtDNA (Fig. 7b). As expected, MELAS patient cells showed higher scores for tRNA genes, and cell-type-specific analysis further revealed tRNA scores and pathogenic mt.3243A>G variant heteroplasmy to be highly correlated (Supplementary Fig. 11e, f). Interestingly, distinct distributions of respective scores were observed among all functional groups. While scmtMPM for complex I, complex IV, and mt-rRNA genes were generally elevated, a more bimodal distribution was observed for complex III and V genes (Fig. 7b). Moreover, the contribution to scwMSS scores for complex III and V variants was generally lower (Fig. 7b; lower panel), likely due to the smaller number and shorter length of their genes, as well as their overall relatively lower MLC scores compared to other mtDNA gene groups (Supplementary Fig. 11g, h)42. These findings indicate distinct mutation dynamics, as also previously intimated across cancers, where complex I subunits were enriched for pathogenic mutations, whereas complex V mutations were broadly depleted for all non-synonymous mutations6.

Upon further evaluation of elevated scMMB scores across distinct functional groups, we followed up on the unexpected observation of high scmtMPM and scwMSS scores in complex I for the pediatric donor (H05) and two MELAS patients (M60 and M80). Specifically, individual gene-level analysis highlighted aberrant increases in scMMB for MT-ND4L in H05 and MT-ND3 in both M60 and M80 (Fig. 7b, c). To identify the underlying mtDNA variants, we ranked all detected variants by their pseudo-bulk heteroplasmy. Notably, in each instance, only a single variant exceeded 1% VAF: mt.10599G>A in MT-ND4 in the H05 (mean VAF = 8%), mt.10270T>C in MT-ND3 in M60 (mean VAF = 25%), and mt.10398A>G in MT-ND3 in M80 (mean VAF = 5.8%), respectively (Fig. 7c, d; Supplementary Table 1). The mt.10599G>A variant in H05 showed a broad and uniform representation among immune cells, including T, NK cells, and monocytes, marking 22% of the entire population of cells (Fig. 7d). Given the low likelihood of such a pronounced clonal expansion at the donor’s young age, this suggests the mt.10599G>A variant to have already been present at higher VAF in the zygote, upon which it further distributed to multiple multipotent hematopoietic stem cells during development as recently suggested by the whole genome-sequencing based evaluation of cellular phylogenies of hematopoietic colonies3,57. Moreover, we observed pronounced enrichments of cells with high mt.10559G>A VAFs in both CD4+ and CD8+ memory T cells, indicating the mtDNA variant tracks with clonal expansion and memory formation (Fig. 7e). In contrast, the mt.10270T>C variant in M60 exhibited a narrow and exclusive distribution marking 76% of CD8+ T effector memory cells (Fig. 7d). In concordance with a very high MLC score of 0.969, SIFT and AlphaMissense predict mt.10270T>C to be deleterious, which likely explains why this variant has not been reported in ClinVar due to a strong purifying selection at the population level22. However, the strong enrichment of the mt.10270T>C mutant allele in memory CD8+ T cells suggests substantial clonal expansion and lineage-specific tolerance toward a pathogenic variant (Fig. 7e). This specific population also notably deviated from the regression line of the scmtMPM and scwMSS (Supplementary Fig. 10b; red arrow), indicating that while both scores are generally correlated, scwMSS integrates the pathogenicity of variants that are not fully reflected by the overall mutational burden alone. The mt.10398A>G variant in M80 is present in various cell types, with 25% of CD8+ naive T cells marked (Fig. 7d), again suggesting the occurrence of likely clonal events associated with individual mtDNA mutations. This variant has been reported to be pleiotropic and associated with distinct disease phenotypes in MitoMAP, such as Leigh syndrome, breast cancer, metabolic syndrome, and Parkinson’s disease58. While this variant appeared to have medium-to-low VAFs in most immune cell populations, we observed enrichment for cells with high mt.10398A>G VAFs in monocytes and CD8+ T cell effector/memory populations (Fig. 7e), again suggesting a potential lineage-specific preference for individual mtDNA mutations. In summary, the evaluation of cellular scMMB scores readily identifies functionally relevant complex- and gene-specific variants, providing insight into their distribution, potential pathogenicity, and implications for cell lineage and mitochondrial function, and motivating further functional investigation.

Methods
Our research complies with all relevant ethical regulations, including ethics approval and guidance from Charité—Universitätsmedizin Berlin, the Berlin Institute of Health at Charité, and the Max Delbrück Center for Molecular Medicine (MDC), Berlin Institute for Medical Systems Biology (BIMSB), Berlin, Germany.

Healthy donor and patient samples
Primary human peripheral blood samples were collected after informed consent and approval was granted by the local ethics commission at Charité—Universitätsmedizin Berlin, Germany (EA1/192/20, EA2/144/15, EA2/284/25). PBMCs were isolated using Ficoll Paque Plus (Sigma Aldrich, GE17-1440) by density gradient centrifugation using SepMate tubes (StemCell Technologies) and were stored in vapor-phase liquid nitrogen after cryopreservation with FBS supplemented with 10% dimethyl sulfoxide until the day of the single-cell sequencing experiment.

Cell culture
POLGD274A knock-in cell lines, with (Knock-in clone A2, KIA2) or without (Knock-in clone 36, KI36) the introduction of a tetracycline-inducible mitochondrial restriction exonuclease mitoEagI, were previously described23. The knock-in line KI36 was derived from the parental control cell line (CTRL), with the KIA2 line being derived from more distantly related HEK293T cells. Briefly, control and KI36 HEK293 cell lines were cultured in high-glucose Dulbecco’s Modified Eagle Medium (DMEM; Gibco, 41966052) supplemented with 10% [v/v] FBS Standard, 50 µg/L Uridine, and 1 U/ml Pen-Strep. KIA2 HEK293T cell lines were supplemented additionally with 50 µg/ml Hygromycin B, and 15 µg/ml Blasticidin S Hydrochlorid without the induction of mitoEagI. When the cell culture reached 90% confluency, cells were detached with Trypsin-EDTA (Gibco, 25200056) and split at 1:10 or cryopreserved. Cells were cultured for up to 2 weeks and underwent up to three passages before single-cell sequencing to limit the accumulation of additional mtDNA mutations.

Galactose treatment
Control and POLGD274A knock-in cell lines were seeded in high-glucose DMEM and incubated at 37°C in a 5% CO₂ incubator overnight to allow cell attachment. The next day, cells were rinsed with PBS, and the medium was replaced with glucose-free DMEM (Gibco, 11966025) supplemented with 1 mM sodium pyruvate (Gibco, 11360070) and with either 10 mM glucose (3 days) or 10 mM galactose (1 day or 3 days) (ThermoFisher, J15995.A1), as indicated. On the day of single-cell sequencing, cells were washed with PBS to remove dead cells, detached with trypsin, and then washed with PBS. Each cell line–treatment combination was labeled with dedicated TotalSeq-B hashtag antibodies, washed, and pooled at equal ratios. The pooled cells were subsequently processed using the mtscATAC-seq protocol as described below.

Mitochondrial enzyme activity measurements
The activities of the mitochondrial enzymes NADH:CoQ oxidoreductase (complex I, CI), cytochrome c oxidase (complex IV, CIV), and citrate synthase (CS) were measured as previously described69. Briefly, cells were homogenized by sonication. Complex I activity was measured as the rotenone-sensitive rate of NADH oxidation in a medium containing 50 mM KCl, 1 mM EDTA, 10 mM Tris-HCl (pH 7.4), 2 mM KCN, 100 µM CoQ1, and 150 µM NADH. Total NADH oxidation was recorded for 2 min before adding 20 µM rotenone to determine the rotenone-insensitive rate. Complex I activity was defined as the rate difference between total and rotenone-insensitive NADH oxidation. Complex IV activity was measured in 0.1 M potassium phosphate buffer (pH 7.4) containing 0.02% lauryl maltoside, and 200 µM ferrocytochrome c (reduced cytochrome c), and ferrocytochrome c oxidation was monitored at 510/535 nm. All measurements were performed at 30°C using a dual-wavelength spectrophotometer (Aminco DW-2000, SLM Instruments, Rochester, NY, USA).

Cell hashing and FACS sorting
Cryopreserved cell lines were thawed in RPMI with 10% FBS by serial dilution. TotalSeqTM-A anti-human hashtag antibodies were purchased from BioLegend. Two million cells per cell line were resuspended in cell staining buffer (BioLegend, 420201) with Human TruStain FcX (BioLegend, 422301) for 10 min on ice. Cell lines were then stained with 0.25 µg of individual hashtag antibodies (see Supplementary Table 2 for barcode assignment) in the cell staining buffer. After staining, cells were washed three times with PBS/BSA. All three HEK293 cell lines were pooled at approximately equal ratios. For live/dead cell discrimination, Sytox Blue (ThermoFisher, S34857) was used at a 1:1000 dilution in PBS/BSA/EDTA. Cell sorting was conducted using the BD Bioscience FACSARIA III cell sorter with a 100-µm nozzle at the MDC/BIH Genomics Platform.

mtscATAC-seq with cell hashing
mtscATAC-seq libraries were generated as previously described and according to the manufacturer’s instructions (CG000209-RevF) with the Chromium Controller and Chromium Next GEM Single Cell ATAC Reagent Kits v1.1 (PN-1000176) with minor modification and sample multiplexing with TotalSeqTM-A hashtag antibodies. Briefly, sorted cells were fixed in 1% formaldehyde (FA; ThermoFisher) in PBS for 10 min at room temperature, then quenched with glycine solution to a final concentration of 0.125 M before washing twice with PBS/BSA by centrifugation at 400g for 5 min at 4°C. Cells were subsequently permeabilized with lysis buffer (10 mM Tris-HCL pH 7.4, 10 mM NaCl, 3 mM MgCl2, 0.1% NP40, 1% BSA) for 3 min (primary cells) or 5 min (cell lines) on ice, followed by adding 1 ml of chilled wash buffer (10 mM Tris-HCL pH 7.4, 10 mM NaCl, 3 mM MgCl2, 1% BSA) and inversion before centrifugation at 500g, 5 min, 4°C. The supernatant was discarded, cells were diluted in 1x Diluted Nuclei buffer (10x Genomics), counted using Trypan Blue and a Countess III FL Automated Cell Counter, and processed for Tn5-based tagmentation. After tagmentation, the cells were loaded onto a Chromium controller to generate single-cell Gel Bead-In-Emulsions (GEMs), followed by linear PCR as described in the protocol, using a C1000 Touch Thermal cycler with a 96-Deep Well Reaction Module (Bio-Rad). To enable the amplification of the Totalseq-A hashtag antibody, 0.5 µl of 1 µM bridge oligo A (BOA) was added to the barcoding mix. Subsequent ATAC and cell-hashing library generation were conducted as described in the user manual and as also outlined online at https://cite-seq.com/asapseq/, respectively.

Sequencing and demultiplexing
All libraries were sequenced using the Illumina NovaSeq 6000 sequencing platforms at the Genomic Core Facility of BIH/MDC (Facility ID = 1565, The Core Marketplace: MDC & BIH Technology Platform Genomics). mtscATAC-seq libraries were sequenced with a 2×100 paired-end read configuration. Raw-sequencing data were demultiplexed using CellRanger-ATAC ‘mkfastq‘.

scATAC-seq analyses
Raw fastq files/sequencing reads were aligned to the mtDNA blacklist modified hg38 reference genome using CellRanger-ATAC count (v2.1.0). Peak-cell matrices and fragments files were further preprocessed using the Seurat v.4 and Signac v.1.10 R packages. Doublet detection was performed using AMULET v.1.1 with default parameters (https://github.com/UcarLab/AMULET)70. For quality control, cells meeting all the following criteria were retained for downstream analysis: >1.000 but <25.000 unique nuclear fragments and a TSS enrichment score >3. Normalization and dimensionality reduction were performed using latent semantic indexing (LSI). The cell hashtag abundances were quantified using ASAP_to_kite.py (https://github.com/caleblareau/asap_to_kite), followed by the kallisto, bustools, and kite pipelines. Cell hashtag counts were then normalized and scaled using the centered log-ratio (CLR) transformation and the ScaleData function implemented in Seurat. For PBMC, cell-type annotations and UMAP coordinates were transferred and mapped from the Azimuth CITE-seq reference dataset using public 10x Genomics Multiome (RNA- and ATAC-seq) PBMC data as a cross-modality bridge via the Seurat Dictionary Learning bridge integration (v4.0.3) method71. The differentially accessible gene (DAG) analysis of galactose-treated cells was performed on the gene score matrix calculated by the scATAC-seq analysis package ArchR72 (v1.0.3) using the using the build-in ‘getMarkerFeatures‘ with the default Wilcoxon pairwise test method.

mtDNA variants calling
Mitochondrial DNA variants were identified from the ATAC-seq modality using mgatk (v0.7.0) (https://github.com/caleblareau/mgatk). Cells with mtDNA coverage <5 were excluded. mtDNA variant allele frequencies were determined as described in ‘variant_calling‘ (https://github.com/caleblareau/mtscATACpaper_reproducibility/), returning a SummarizedExperiment object, containing two assay matrices, ‘allele_frequency‘ and ‘coverage‘. High-quality mtDNA variants were identified by filtering variants with a strand concordance >0.65, mean-variance ratio >0.01 in the PBMC dataset, mean-variance ratio >0.0031 in the HEK293T dataset, and n_conf_cells_detected ≥ 1, unless indicated otherwise. In addition, four mtDNA variants that are known to be sensitive to sequencing artifacts due to their homopolymer sequence context, including mt.301A>C, mt.310T>C, mt.316G>C, and mt.3302A>C, were excluded from downstream analysis59. For the “relaxed” filtering, we additionally examined mtDNA variants with a strand concordance between 0.45 and 0.65. Notably, these variants retained the expected replication-error mutational signature, as demonstrated in Supplementary Fig. 1f. However, all other analyses presented throughout the manuscript exclusively utilized variants identified via the stricter filtering criteria. The list of pathogenic mtDNA variants was obtained from the MITOMAP database (https://www.mitomap.org/foswiki/bin/view/MITOMAP/ConfirmedMutations). Only variants annotated as “Confirmed (Cfrm)”, as “Confirmed–Pathogenic (Cfrm[P])”, or “Confirmed–Likely Pathogenic (Cfrm[LP])” were retained. The resulting list of confirmed pathogenic variants used for annotation is provided in Supplementary Data 4.

Sequencing downsampling analysis
Position-sorted BAM files generated from the ‘cellranger-atac’ pipeline were processed with samtools (v1.22). Reads aligning to the mitochondrial chromosome (chrM) were extracted, and cell barcode lists obtained from hashtag-based demultiplexing were used to separate individual cell line groups (CTRL, KI36, and KIA2). To normalize sequencing depth across samples, BAM files were randomly downsampled using samtools view -s (seed=42) to fixed read depths of 10M, 20M, 30M, 40M, 50M, 100M, 150M, and 200M reads. The resulting downsampled BAM files were subsequently processed with mgatk (v0.7.0) for mtDNA variant calling. To visualize depth-dependent changes in variant composition, we normalized the pseudobulk VAFs of each mtDNA variant to the sum heteroplasmy of all variants at each downsampled depth (relative heteroplasmy). For the ΔVAF density plots, we first filtered mtDNA variants detected in both the 50M and 200M KIA2 datasets, and computed per-cell, per-position ΔVAF as VAF200M – VAF50M, excluding variants with zero VAF at both depths in that cell. The resulting ΔVAF values were visualized using kernel density estimation to summarize the distribution of single-cell VAF changes between 50M and 200M total mtDNA reads.

Single-cell mitochondrial DNA mutational burden (scMMB)
We adapted metrics previously used to assess mtDNA mutational burden[6] and local constraint22,41,42. These were modified into single-cell mtDNA MPM (scmtMPM) and single-cell Heteroplasmy-Weighted MSS (scwMSS). The mtDNA mutational burden in single cells was calculated at the levels of individual genes, respiratory chain complexes (I, III, IV, and V), and the whole mtDNA genome, as previously conducted in whole-exome sequencing data6. Mutation counts were calculated by first multiplying the ‘allele_frequency‘ and ‘coverage‘ matrices, then summing across the respective mitochondrial genomic regions. The total genomic length sequenced was calculated by summing up the coverage matrix of the respective mitochondrial genomic regions of individual cells. Note that the scmtMPM score is normalized to mtDNA-aligned coverage (i.e., the total number of bases sequenced within the defined mtDNA region per cell). mtDNA mutation per million base pairs (MB) was calculated as the total somatic mtDNA mutation count divided by the total genomic length sequenced (in MB) for the respective genomic regions of a single cell. The weighted mitochondrial local constraint (MLC) score sum (MSS)22 was calculated by multiplying the position-wise MLC score by its single-cell heteroplasmy and summing the weighted scores for the respective genomic regions within a single cell.

Comparison of variant calling from mtDNA- and mtRNA-derived sequencing reads
To assess the utility of popular RNA transcriptomic workflows for mtRNA variant calling, we quantified VAFs from both modalities from a 10x Multiome/DOGMA-seq library of PBMCs from the same healthy pediatric donor11 using a cross-calling and filtering workflow (Supplementary Fig. 12a). First, variants were primarily called from mtDNA using mgatk to exclude artifactual variants that have previously been shown to be enriched in RNA-seq data. Second, variants identified from mtDNA were cross-called in mtRNA reads. To account for uneven mtRNA transcript coverage, we further filtered the cell-by-position mtDNA/mtRNA coverage matrices and retained only positions with at least five reads per cell. We then filtered and re-quantified mutation counts and VAFs matrices with the filtered coverage matrices and compared the number of confidently detected variants in both mtDNA and mtRNA (Supplementary Fig. 7c). Modality weight was calculated as the number of detectable mtDNA variants divided by the sum of detectable variants in mtDNA and mtRNA (plus 0.01 for computational reasons).

Statistics and reproducibility
No statistical method was used to predetermine the sample size. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessments. Two-sided Wilcoxon rank-sum tests were performed, and P-values were adjusted for multiple comparisons using the Holm–Bonferroni procedure, with significance defined at an FDR of 5%. All statistical analyses were conducted in R, and custom code used to replicate the analyses is available as part of the code availability.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is an interesting paper which presents a very thorough bioinformatics analysis of cell lines with an inducible mutation in\nPOLG (KI) compared to control lines, followed by an analysis of human patient blood cells from individuals with know\npathogenic mtDNA mutations. The main novel finding is a much greater number of mtSNVs in the POLG KI cells than\nanticipated, and an increase in mtSNV rate in the human blood cells. This is supplemented by a detailed annotation of the\nmtSNVs which show a pattern consistent with selection against functional variants (which is expected and confirmatory).\nThe technical approach is cutting edge and the analysis thorough, with a very sensible novel development of variant\nconstraint analysis which encompasses heteroplasmy.\n- It is not always explicit whether the analysis is based on data at single cell resolution or in pseudobulk. It would help to\nshow the cell-level mtDNA coverage. I infer this to be very in controls low because the 2-3 fold increase described in POLG\nKI mutant cells was <30x, so does this mean the control coverage was <<10x? The current draft presents normalised\ncoverage which conceals this raw data. This is very important because low coverage will drastically limit the ability to detect\nand measure heteroplasmy levels accurately – introducing a systematic bias between the cell lines which will not be\naddressed in a pseudobulk analysis. I am afraid this is potentially a fundamental flaw.\n- The POLG KI line is inducible allowing the authors to turn on and turn off the mutator phenotype. This is very elegant but\nnot capitalised on. Why didn’t the authors compare the same cell line before and after a burst of mutagenesis? This would\ncontrol for many of the confounders they struggle with (eg inherited variants – and critically, the method used to generate the\ncell line itself which could induced mutations, or lead to rapid segregation of pre-existing mutations through bottleneck\neffects). I am afraid this could also be another fundamental flaw influencing the main result.\n- The hash-tagging method is not without difficulties but the authors do not show that their method reliably separates the\ndifferent cell lines. This is critically important. We need to be confident that the authors are able to reliably resolve single cell\nsequences to measure the mutation burden in single cells.\n- Orthogonal validation. The methods they are using are novel and it would be reassuring to see some orthogonal validation\neg of the differences in mtDNA content (related to depth), and the presence of linearised mtDNA at high levels. Genuine bulk\nanalysis of the mtDNA, deletion and transcriptomes would provide this.\n- It would be extremely reassuring to see the effects of growing the cell lines in galactose, which forces oxidative\nmetabolism. I would anticipate that the selection signal would be much stronger – validating their conclusions\nexperimentally, rather than by inference from single observations in each cell line.\n- What is the provenance of the cell lines? Eg number of passages for the primary and POLG KI lines?\n- I am unclear whether their pipeline excludes the confounding effects of mtDNA sequences embedded in the nucleus\n(important for the low-heteroplasmy variants)\n- Definition of de novo variants. As the authors note, there is a threshold for detection of mtSNVs linked to mtDNA coverage. It\nis therefore impossible to be sure whether any variant is strictly ‘de novo’ or whether it was present in a small number of\ncells at low heteroplasmy levels. The authors must temper their conclusions/results with this in mind.\n- Related to this – the ability to detect low level heteroplasmies is dependent on the sequencing depth. Their own data\nshows site-specific differences in coverage. Did the authors account for this in their analysis? Particularly when comparing\ndifferent mutation classes (eg could it be that pathogenic sites tend to have low coverage, thus limiting their potential for\ndetection). It would be useful for the field to show whether this is an important confounder of not.\n- The authors have previously used mtDNA to lineage trace, but they do not do this in the current paper, despite their\nassertion that there are sub clones within their cell lines. I think this needs greater exploration before they reach the\nconclusions that they do.\n- Finally, the manuscript refers to phenotypes on several occasions, but their current analysis does not systematically study\nphenotype. Eg ‘these findings clearly suggest POLGD274A phenotypes to be driven by the cumulative mutational burden as\ncompared to the impact of individual or a few deleterious variants as observed in classical mitochondriopathies.’ I am afraid\nthis conclusion cannot be drawn from the data presented which does not include these phenotypes and certainly doesn’t\nshow a causal role.", "reviewer_2": "In this manuscript by Hsieh et al, the authors leverage a technique mitochondrial single-cell ATAC-seq, which they\npreviously established to investigate heteroplasmy levels of specific variants from patients with mtDNA disease, (Nature\nGenetics 2022) to measure mtDNA mutational loads and somatic mosaicism. To bench-mark their technique they use HEK-\n293 cell lines with a knock-in mutation in POLG which causes a proof-reading deficiency of the mtDNA polymerase and\nelevated mtDNA mutation rates. Interestingly in these lines, the variant allele frequencies (VAFs) of individual cells did not\nreach classical thresholds for metabolic dysfunction. However, when complex I and complex IV activity was measured in\nbulk cell cultures, there was metabolic impairment, suggesting multiple mutations can act synergistically to affect\nmitochondrial function.\nThe authors then introduce two metrics: single-cell mtDNA mutations per million base pairs (scmtMPM) and heteroplasmy-\nweighted mitochondrial local constraint scores (scwMSS). These useful measures allow comparative analysis of samples\nfrom subjects of different ages and disease groups which is very useful for the research community. They then apply their\nbenchmarked technique and measures to test datasets including control blood samples from young (5 year old) and older\n(47 year old) controls and 4 patients with the m.3243A>G mutation (aged 29-80). They show that the somatic mtDNA\nmutation load per cell increases with age in control samples and suggest that the m.3243A>G mutation may affect the\naccumulation of somatic mtDNA mutations.\nSpecific comments:\n1. The authors have established a protocol to quantify mtDNA mutation load at the single cell level at an unprecedented\nresolution. Previous single cell work on somatic mtDNA mutations has been at the level of 100s of cells, this was 1000s. The\ndata is very high resolution and supports previous studies which have shown similar mutational spectra and dynamics of\nage-related somatic mtDNA mutation accumulation.\n2. This work is of high technical quality and provides a robust framework for assessing mtDNA load/constraints/evolution in\nother tissue/disease subgroups.\n3. I think the title could be clearer, as although the technique is technically multi-omic, other than one experiment towards the\nend where the nuclear ATACSeq data is used to subset immune-cells, this paper is really focussed on mtDNA sequence\nanalysis. There is very little discussion of the impact (if there is any) of cell type on somatic mtDNA mutagenesis. The\nmajority of the analysis groups all of the different cell types together. I would suggest either fully describing all of the data at\nthe level of cell types, or just using the name of the technique in the title (mitochondrial single-cell ATAC-seq) which will\nshow the focus of the paper.\n4. There appears to be three distinct populations of cells in the KI36 cell line, but the data is analysed as one population. Are\nthere specific clonal variants that are driving the clustering in this line?\n5. Can you be more explicit when describing the confirmed pathogenic mutations on Mitomap? There are a number of\ndifferent tables on there and I’m guessing you used the ‘confirmed’ table. It would be good to state this to avoid confusion.\n6. For the scmtMPM score, from what I understand, this is normalised to coverage? Can you be explicit that this is the total of\nthe mtDNA aligned reads and not the overall sequencing depth of the ATAC-Seq?\n7. I don’t really understand how the conclusion was drawn that the 3243 mutation impedes somatic mtDNA mutation\naccumulation? In figures 5d and 5e both of these patients are relatively young (29 and 35) and the total scmtMPM\n(-3243A>G) looks to be slightly lower than H47 (closest age-match). However, M60 and M80 are much higher supporting the\naccumulation of somatic mtDNA mutations with age. What is the evidence that the 3243 variant impedes the accumulation of\nsomatic variants? What happens if you draw similar plots for M60 and M80? I may be misunderstanding the figure, is there\nan alternative way of presenting this data?", "reviewer_4": "The study by Yu-Hsin Hsiehand colleagues presents a study in mitochondrial genomics through the development and\napplication of novel single-cell metrics to quantify mitochondrial DNA (mtDNA) mutational burden and heteroplasmy\ndynamics. Utilizing mitochondrial single-cell ATAC-seq (mtscATAC-seq), the authors introduce two key metrics—single-cell\nmtDNA mutations per million base pairs (scmtMPM) and heteroplasmy-weighted mitochondrial local constraint scores\n(scwMSS)—to capture genome-wide mtDNA variation at single-cell resolution. The study benchmarks these metrics in\nPOLGD274A knock-in cell lines and applies them to primary human immune cells from both healthy donors and patients\nwith mitochondrial encephalomyopathy (MELAS). While the work addresses an important challenge in mitochondrial\ngenetics, offering potentially valuable insights into mitochondrial disease mechanisms, several critical issues and concerns\nwarrant further attention.\nMajor Comments:\n1. The authors reported 455 variants in control cells versus 2484 and 4651 variants in the KI36 and KIA2 POLGD274A\nknock-in cells, respectively. Given the substantially higher sequencing depth in the knock-in lines, the increased variant\ncounts may largely reflect depth-related bias. Additionally, the authors should also clarify how uneven mtDNA coverage\nacross the mitochondrial genome was accounted for during variant calling and downstream analyses.\n2. The two metrics (scmtMPM and scwMSS) represent key contributions of this study. However, biases in variant detection\nand variability in sequencing coverage could substantially influence their accuracy. The reported median mtDNA\nsequencing depth, ranging from approximately 10x to 40x, is very low for confident variant calling, particularly for low-\nfrequency variants. Although the authors have implemented measures to minimize these biases, further validation, such as\ntargeted deep sequencing, is needed to confirm the accuracy of variant calls and to robustly support the conclusions drawn\nfrom these metrics.\n3. The analysis assumes that variants unique to POLGD274A lines are de novo, but clonal expansions or bottlenecks could\nalso influence variant distributions. Could the authors elaborate on how these factors were considered or controlled for in\ntheir analyses?\n4. The authors incorporate a mitochondrial genome constraint model based on population genetics data from gnomAD.\nHave the authors considered the impact of population diversity, mitochondrial haplotypes, as well as the cell-type-specific or\ncontext-dependent functional constraints?\n5. Have the authors explored integrating other variant annotations, like predicted pathogenicity scores or conservation, to\nhelp refine the MLC-based scoring?\n6. Given that mtscATAC-seq simultaneously captures mtDNA variants and chromatin accessibility, have the authors\ninvestigated potential interactions between mtDNA mutational burden and nuclear chromatin states? Such as, whether/how\nnuclear-mitochondrial interactions influencing variant dynamics? whether cells with higher mutational burden exhibit altered\nchromatin accessibility; Additionally, single-cell Multiome technique, which profile both transcriptomes and chromatin\naccessibility in individual cells, could provide an even more powerful platform to investigate the functional consequences of\nmtDNA mutations. Have the authors considered applying or integrating such approaches?\n7. It would strengthen the impact of the study if the authors included functional assays to see how mtDNA variants actually\naffect cells. Further linking the scmtMPM and scwMSS metrics to cellular phenotypes would add important biological context\nto the current study.\n9. Could the authors expand on the translational implications of their findings. In particular, discussing how feasible it might\nbe to apply the scmtMPM and scwMSS metrics in clinical or therapeutic settings, along with any potential challenges, would\nadd valuable perspective.\n9. The manuscript currently lacks enough detail about the statistical analyses used. The authors should include clear\ndescriptions of the statistical methods applied throughout the study to better support their findings and improve transparency.\nOther Comments:\nIt would be helpful if the authors provided a summary table listing all variants detected in this study, including average\nsequencing depth, the number of cells carrying each mutation within each group, and other relevant information."}, "author_response": "We thank the Reviewers for the very thorough review of our work, which helped to significantly\nimprove our manuscript. We provide a detailed point-by-point response to all comments in blue and\nhave modified the manuscript based on the valuable suggestions of the Reviewers. Corresponding\nchanges in the manuscript are also highlighted in blue.\nReviewer #1 (Remarks to the Author):\nThis is an interesting paper which presents a very thorough bioinformatics analysis of cell lines with\nan inducible mutation in POLG (KI) compared to control lines, followed by an analysis of human\npatient blood cells from individuals with know pathogenic mtDNA mutations. The main novel finding\nis a much greater number of mtSNVs in the POLG KI cells than anticipated, and an increase in mtSNV\nrate in the human blood cells. This is supplemented by a detailed annotation of the mtSNVs which\nshow a pattern consistent with selection against functional variants (which is expected and\nconfirmatory). The technical approach is cutting edge and the analysis thorough, with a very sensible\nnovel development of variant constraint analysis which encompasses heteroplasmy.\nWe thank the Reviewer for their praise and for acknowledging the advances presented in our\nmanuscript.\n1. It is not always explicit whether the analysis is based on data at single cell resolution or in\npseudobulk. It would help to show the cell-level mtDNA coverage. I infer this to be very in controls low\nbecause the 2-3 fold increase described in POLG KI mutant cells was <30x, so does this mean the\ncontrol coverage was <<10x? The current draft presents normalised coverage which conceals this\nraw data. This is very important because low coverage will drastically limit the ability to detect and\nmeasure heteroplasmy levels accurately – introducing a systematic bias between the cell lines which\nwill not be addressed in a pseudobulk analysis. I am afraid this is potentially a fundamental flaw.\nWe thank the Reviewer for these important comments. We agree that low mtDNA coverage can limit\nthe ability to detect and accurately measure heteroplasmy. To investigate this, we generated 8x-fold\nmore deeply sequenced data for our HEK293 lines, which improved the raw reads from about 23,000\nto 190,000 per cell (10x Genomics recommends up to 25,000 reads) and yielded a substantially\nincreased per-cell mtDNA depth (median: CTRL = 55.6×, KI36 = 154.1×, KIA2 = 255.2×). Importantly,\nthe increased mtDNA depth improved variant detection sensitivity and robustness; however, the\nrelative differences between control and POLGD274A knock-in lines remained consistent with our\nprevious results. Likewise, the scmtMMB metrics remained stable across varying sequencing depths\n(new Extended Data Fig. 7), further supporting the robustness of our approach.\nAs such, we updated our HEK293 results with the more deeply sequenced data in our revised\nmanuscript (Figs. 1–4) and further clarified throughout the text whether analyses are performed at\nthe single-cell or pseudobulk level. In addition, we would like to clarify that per-cell mtDNA coverage\ndistributions are shown in Fig. 1c (projected onto the UMAP), the total coverage being quite uniform\n(Fig. 1e) and Extended Data Fig. 1d (as a violin plot), to which we now also directly refer in the main\ntext and figure legends. We have provided these figure panels also below as Reviewer Fig. 1 for the\nReviewer’s more direct consideration.\nReviewer Fig. 1 – Raw mitochondrial DNA (mtDNA) depth in single cells. UMAP projection color-scaled by\nthe mtDNA depth (left), the violin plot of raw mitochondrial DNA depth (middle), and the mtDNA position-wise\npseudobulk coverage (right).\nIn addition, we revised our discussion to emphasize that a minimal, and importantly comparable,\nmtDNA sequencing depth is required for meaningful biological comparisons. We further point out that\nbiological differences in mtDNA copy number across cell types and analytical methods may need to\nbe accounted for when interpreting such results.\n2. The POLG KI line is inducible allowing the authors to turn on and turn off the mutator phenotype.\nThis is very elegant but not capitalised on. Why didn’t the authors compare the same cell line before\nand after a burst of mutagenesis? This would control for many of the confounders they struggle with\n(eg inherited variants – and critically, the method used to generate the cell line itself which could\ninduced mutations, or lead to rapid segregation of pre-existing mutations through bottleneck effects).\nI am afraid this could also be another fundamental flaw influencing the main result.\nWe thank the Reviewer for this thoughtful comment. We would like to clarify that the inducible element\nin the POLGD274A knock-in KIA2 line refers to the tetracycline-inducible mitoEagI restriction enzyme,\nwhich cleaves mtDNA but does not affect the proofreading deficiency due to the POLGD274A mutation\nitself. Thus, mitoEagI induction does not generate a burst of mutagenesis but instead results in mtDNA\ncleavage. However, this is a feature that we did not explicitly utilize in our study. Instead, the POLG\nanalysis was conducted with the primary intention of quantifying differences in the mutational burden\narising from the POLG exonuclease deficiency (as opposed to the kinetics of the system). As such,\nwe used the additional KIA2 line without induction as an additional knock-in sample to benchmark our\nmutational burden analysis.\nTo evaluate selection dynamics of mutation under a metabolic bottleneck, we have now included new\ndata where the cell lines were cultured in the presence of glucose or galactose. Notably, the KIA2 line\nwith the higher mutational burden did not tolerate galactose treatment, prohibiting further evaluation.\nFor the KI36 line, about 75% of cells quickly succumbed in the presence of galactose, suggesting a\nsignificant bottleneck. Notably, instead of pronounced shifts in heteroplasmy, we observed a\nstatistically significant increase in mtDNA copy number, that we interpret as the primary mode of\ncompensation in this experimental setting. We present these results in new Fig. 5 and Extended\nData Figs. 8-9.\nHowever, we believe that different concentrations and/or longer treatments may bring about distinct\nmutation kinetics. We also note that relatively few confirmed pathogenic mtDNA mutations were\npresent in the lines to begin with, which likely would have exerted more pronounced selective\npressures. Also, the mutational burden was relatively uniform across the KI36 line, potentially\nindicating that the effects of galactose treatment may have also uniformly affected cells as opposed\nto selecting for a minor and more tolerant (clonal) population. We believe additional investigations to\nbe necessary and a fruitful avenue to further dissect POLG mutator biology, but we would respectfully\nconsider such experiments to be more suitable for future investigation.\nWe also acknowledge that some residual uncertainty may remain regarding pre-existing mutations,\nin particular those of low heteroplasmy, and present in only a small proportion of cells. However, such\nmutations will always be challenging to reliably capture, given detection limits and or sampling bias\nof all sequencing assays. As such, we would not consider our result to be fundamentally flawed as\nopposed to being subject to typical experimental and analytical restrictions of this type of data.\nAnalogous to single-cell transcriptomic analysis, the capture of lowly expressed transcripts or small\ncell populations can be biologically meaningful, but care should be taken with respect to biological\ninterpretation, as we now more explicitly point out limitations in our revised discussion.\nGiven our primary goal of benchmarking single-cell assessments of mtDNA mutational burden and\nestablishing appropriate quantitative metrics, we believe that our more deeply sequenced HEK\ndataset, the new galactose experiment, together with the additional revisions to the manuscript,\nconstitute a significant and robust methodological advance for the field.\n3. The hash-tagging method is not without difficulties but the authors do not show that their method\nreliably separates the different cell lines. This is critically important. We need to be confident that the\nauthors are able to reliably resolve single cell sequences to measure the mutation burden in single\ncells.\nWe thank the Reviewer for this comment. Hashtag-based multiplexing is a widely used and well-\nvalidated approach for single-cell demultiplexing1,2, and we followed these established protocols and\nanalytical workflows. In our dataset, the three HEK293 cell lines showed clear and distinct hashtag\nsignal separation, as illustrated in new Extended Data Fig. 1c, demonstrating that the hashtagging\nmethod reliably resolved the individual cell lines for confident downstream single-cell mutation burden\nanalyses.\nIn our revised manuscript, we again applied hashtag-based multiplexing in a more complex setting.\nIn the galactose treatment experiment, we use six hashtag antibodies to distinguish the cell line-\ntreament combinations (CTRL and KI36 cultured in glucose for 3 days, or galactose for 1 and 3 days).\nAs shown in new Extended Data Fig. 8b-c, the six hashtags show well-separated signal distributions\nwith minimal overlap, supporting the robust demultiplexing of all conditions and validating the\nperformance of the methodology as previously benchmarked.\n4. Orthogonal validation. The methods they are using are novel and it would be reassuring to see\nsome orthogonal validation eg of the differences in mtDNA content (related to depth), and the\npresence of linearised mtDNA at high levels. Genuine bulk analysis of the mtDNA, deletion and\ntranscriptomes would provide this.\nWe thank the Reviewer for this comment. Regarding the linearized mtDNA, we have previously\ndemonstrated its detection via an orthogonal Southern blot, and next-generation sequencing of linker-\nligated bulk mtDNA isolated from immunoprecipitated mitochondria3.\nRegarding the mtscATAC-seq and variants assay thereof, we emphasize that these have been\nextensively benchmarked and applied2,4–8. Modified versions have further been adapted by members\nof the mitochondrial genetics’ community, such as for MitoPerturb-seq9. Earlier work demonstrated\nthat single-cell heteroplasmies were further well correlated across different single-cell genomic\nassays, including scRNA-seq and scATAC-seq4,5. In that work, we also demonstrated that pseudobulk\nheteroplasmies called from single-cell assays are highly correlated with bulk heteroplasmy\nmeasurements from bulk genomic sequencing. Please see, for example, Fig. 2F and Fig. S2E from\nLudwig et al.4, which we pasted below for the Reviewer’s consideration (Reviewer Fig. 2).\n[editorial note: third party material redacted]\nWe further conducted downsampling analysis of the POLGD274A knock-in line data to assess the\ndetection of mtDNA variants as a function of their allelic frequencies at variable mtDNA sequencing\ndepth. As expected, this is very robust for higher allelic frequency variants (new Extended Data Fig.\n7). Importantly, also low-heteroplasmy variants (those <1% VAF) were consistently recovered, and\ntheir pseudobulk VAF estimates remained highly concordant across variable sequencing depths\n(Extended Data Fig. 7f). Moreover, single-cell mtDNA heteroplasmy estimates were highly\nconsistent, with most per-variant ΔVAF values within ±5% and a global ΔVAF distribution tightly\ncentered at zero when comparing a 4-fold difference in sequencing depth (Extended Data Fig. 7g).\nTogether, these analyses suggest that allelic frequency estimates are largely stable once moderate\nsequencing coverage is achieved.\nFurther, we note that all HEK293 lines, including the control, would be expected to have the same\nbackground rate of potential false-positive mtDNA variants. As can be seen in Fig. 2d, the number of\nadditional variants detected in the POLGD274A knock-in line is, however, substantially larger than in\nthe control line. As such, it is extremely unlikely that this population and the substantial increase of\ndetected variants within is confounded by false-positive variant calls.\nWe also note our previously introduced paired-end, strand-aware variant calling workflow, as\nimplemented in the mgatk pipeline, provides higher specificity and more conservative filtering than\nstandard bulk mtDNA workflows5. As such, we would be uncertain to what extent bulk analysis can\nprovide orthogonal validation beyond our previously published evaluations4. In particular,\ntranscriptomic data are subject to various expression levels and, as such, depths of different mtDNA-\nencoded genes, with tRNA transcripts not being typically captured (Extended Data Fig. 12). Further,\nwe previously identified a significant number of false-positive variant calls, which we speculated is\ndue to RNA polymerase transcriptional error (see Reviewer Fig. 1b above, which shows variants\ndetected at the RNA, but not DNA level; some of which being due to RNA editing4).\n5. It would be extremely reassuring to see the effects of growing the cell lines in galactose, which\nforces oxidative metabolism. I would anticipate that the selection signal would be much stronger –\nvalidating their conclusions experimentally, rather than by inference from single observations in each\ncell line.\nWe thank the Reviewer for this comment. As suggested, we cultured the HEK293 control and\nPOLGD274A knock-in lines in the presence of glucose or galactose. Notably, galactose treatment was\nlethal to the KIA2 line, and a large fraction of about 75% of cells from the KI36 line quickly died (new\nExtended Data Fig. 8), suggesting a significant bottleneck. As elaborated above and discussed as\npart of the Reviewer’s prior comment #2 and newly added data in the new Fig. 5 and Extended Data\nFigs. 8-9, we did not observe pronounced shifts in heteroplasmy but instead a statistically significant\nincrease in mtDNA copy number, suggesting the primary mode of compensation in this experimental\nsetting.\nAs discussed above, we believe that varying the concentration or duration of treatment could result\nin different mutation dynamics. As our cell lines contained relatively few confirmed pathogenic mtDNA\nmutations, the associated selective pressures were likely limited. Moreover, the relatively uniform\nmutational burden observed across the KI36 line suggests that galactose treatment may have affected\nall cells similarly, rather than selecting for a more tolerant subpopulation. While we consider further\nstudies in this direction to be promising, we view them as bestsuited for future independent\ninvestigations. However, the significant cell death of the POLGD274A knock-in lines provides strong\nadditional evidence that the elevated mutational burden significantly compromised mitochondrial\nfunction.\n6. What is the provenance of the cell lines? Eg number of passages for the primary and POLG KI\nlines?\nFor the generation of the POLGD274A knock-in lines, single-cell cloning and about two months of\nculturing were required to expand the populations, during which mtDNA mutations accumulated,\nfollowed by cryopreservation of aliquots. Upon thawing, aliquoted cells were cultured for up to three\npassages (less than two weeks) before single-cell sequencing to minimize the accumulation of\nadditional mtDNA mutations in vitro. The passage numbers at the time of sequencing were P18, P19,\nand P16 for CTRL, KI36, and KIA2, respectively. We have now revised our method section to\nincorporate this relevant information.\n7. I am unclear whether their pipeline excludes the confounding effects of mtDNA sequences\nembedded in the nucleus (important for the low-heteroplasmy variants)\nWe thank the Reviewer for this important comment, which we have previously critically evaluated5.\nSince NUMTs are typically not enriched in regions of accessible chromatin, they are not readily\namplified by scATAC-seq. In line with this, in the original mtscATAC-seq study by Lareau et al., we\nhave shown that NUMTs only negligibly contribute to single-cell ATAC-seq-based library preparation,\nexpecting only ~1 fragment per cell. This makes NUMT very unlikely to meaningfully confound\nheteroplasmy estimates. We have pasted the description of our prior evaluation5 directly below for the\nReviewer’s consideration:\n“To estimate the number of accessible NUMT fragments that would be assigned to mtDNA, we\nconsidered two different approaches. First, we used a public GM12878 dataset from 10X Genomics\n(https://www.10xgenomics.com/solutions/single-cell-atac/) that was aligned to the standard hg19\nreference and counted the number of fragments per cell overlapping our NUMT blacklisted regions,\nwhich resulted in mean 1.4 and median 1.0 fragments per cell. Second, we used a compendium of\nDNase accessible peaks from 164 distinct samples from the ENCODE and Roadmap consortia, and\nestimated that these samples contained a mean 22.6 peaks overlapping our NUMT blacklist. Next,\nusing the GM12878 peak set and the same scATAC-seq dataset, we determined that mean 4.1% of\nthe GM12878 DNase peaks were detected over all cells. The product of these two numbers\n(22.6 × 0.041 = 0.93 fragments per cell) provides an alternative estimate for the number of accessible\nchromatin fragments overlapping NUMTs (~1 fragment) that were blacklisted. As our mtscATAC-seq\nassay generates ~5,000–10,000 mtDNA fragments, we conclude that our blacklist approach yields\nnegligible NUMT contamination.”\n8. Definition of de novo variants. As the authors note, there is a threshold for detection of mtSNVs\nlinked to mtDNA coverage. It is therefore impossible to be sure whether any variant is strictly ‘de novo’\nor whether it was present in a small number of cells at low heteroplasmy levels. The authors must\ntemper their conclusions/results with this in mind.\nWe thank the Reviewer for pointing out this important limitation, which we fully agree with. Along these\nlines, there may be sampling biases, in particular for variants that are only present in a very small\npopulation of cells. As such, in the revised manuscript, we have replaced the term “de novo” with\n“POLGD274A-line specific” to accurately reflect the fact that strictly “de novo” variants cannot be\nunequivocally defined due to coverage-dependent detection thresholds and other biases, such as\nsampling. However, we would like to emphasize that, although a strictly “de novo” definition is not\npossible, we would interpret our deeper sequencing and subsampling analyses to indicate that the\nsheer number of additional observed mutation events in the POLGD274A knock-in cell lines suggests\nthat most detected variants are likely to have occurred “de novo” in this genetic background. We have\naccordingly revised our presentation of these results, acknowledging the discussed limitations.\n9. Related to this – the ability to detect low level heteroplasmies is dependent on the sequencing\ndepth. Their own data shows site-specific differences in coverage. Did the authors account for this in\ntheir analysis? Particularly when comparing different mutation classes (eg could it be that pathogenic\nsites tend to have low coverage, thus limiting their potential for detection). It would be useful for the\nfield to show whether this is an important confounder of not.\nWe thank the Reviewer for this thoughtful comment. Our data show near-uniform mtDNA coverage\nacross the mitochondrial genome (Fig. 1e). This observation is consistent with our previous report5,\nwhere we demonstrated near-uniform mtDNA coverage with mtscATAC-seq, with minor residual\nvariation attributed to PCR amplification or Tn5 insertion bias. However, in the POLGD274A knock-in\nlines, we observe a relatively lower coverage across positions 0-6000, reflecting the linearization of\nmtDNA in the POLGD274A knock-in background, as previously reported3. To more directly evaluate\nmtDNA coverage, we plotted the mean coverage stratified by the positions of reported pathogenic\nvariants from MITOMAP, as well as the MLC score interval (Reviewer Fig. 3 below). These analyses\nclearly demonstrate that neither potential pathogenic nor highly constrained sites exhibit any tendency\ntoward lower or more variable coverage. As such, we conclude that mutation classes and mtDNA\nconstraint are unlikely to be subject to systematic confounding due to limitations in coverage or\nmtDNA variant detection.\nReviewer Fig. 3 – Comparable per-base mtDNA mean coverage across pathogenic and constraint\nvariant classes. Violin plots show the per-base mtDNA meancoverage stratified by positions with pathogenic\nvariants annotated in MITOMAP (left) and MLC score intervals [Low (0–0.33), Mid (0.33–0.66), and High (0.66–\n1)] (right) in control and in POLGD274A knock-in HEK293 cell lines.\n10. The authors have previously used mtDNA to lineage trace, but they do not do this in the current\npaper, despite their assertion that there are sub clones within their cell lines. I think this needs greater\nexploration before they reach the conclusions that they do.\nWe thank the Reviewer for this comment. In our manuscript, we have exemplarily showcased\nsubclone-specific mtDNA variants within the KI36 line (Fig. 4e). Using our previously established\npipeline5, we further demonstrate the confident assignment of distinct subclonal groups based on\nshared mtDNA variants (Reviewer Fig. 4). While lineage tracing in POLGD274A knock-in contexts\npresents an exciting direction, we believe it has relatively limited added value to the main focus of the\ncurrent work that is focused on benchmarking and deriving quantitative metrics for single-cell mtDNA\nmutational burden analyses to evaluate the cellular heterogeneity of mtDNA mutational profiles as\ndemonstrated in our manuscript. Moreover, we note that the higher mutation rate in the investigated\nsetting may require additional benchmarking, as mutational events may arise independently\n(equivalent to barcode homoplasy) and thereby confound clonal inference.\nReviewer Figure 4. Subclonal structure of KI36. (a) UMAP of the POLGD274A knock-in HEK293 clone KI36,\ncolored by clonotype ID (1-3), identified using the FindClonotype function in the R-package Signac. Manual\nsubcluster labels (C1–C3) are shown; each clonotype has a one-to-one match with a subcluster. (b) Heatmap\nshowing variant allele frequencies of 188 clonally informative mtDNA variants (rows) in individual KI36 cells\n(columns). Cells are grouped by clonotype (top color bar), highlighting the subclonal structure.\n11. Finally, the manuscript refers to phenotypes on several occasions, but their current analysis does\nnot systematically study phenotype. Eg ‘these findings clearly suggest POLGD274A phenotypes to\nbe driven by the cumulative mutational burden as compared to the impact of individual or a few\ndeleterious variants as observed in classical mitochondriopathies.’ I am afraid this conclusion cannot\nbe drawn from the data presented which does not include these phenotypes and certainly doesn’t\nshow a causal role.\nWe thank the Reviewer for raising this point. We agree that our study does not always directly assess\ncellular phenotypes, as opposed to introducing and benchmarking a broadly applicable single-cell\nframework to quantify mtDNA mutational landscapes. In the original manuscript, we showed altered\ncomplex I and IV activities of the POLGD274A knock-in lines, and our new galactose treatments provide\nadditional evidence of mitochondrial dysfunction attributable to the elevated mtDNA mutational\nburden (new Fig. 5 and Extended Data Figs. 8-9). The pronounced cell death/ growth inhibition of\nthe POLGD274A knock-in lines, elevations in mtDNA copy number in response to galactose, overall\nstable distribution of single-cell VAFs, and alteration in chromatin accessibility profiles of mitochondrial\nstress response-related genes in KI36 under galactose treatment provides additional phenotypical\nevaluation. As such, we have carefully revised the manuscript text to adjust the description of our\nobservations. Our revised version now also acknowledges the need for future investigations and\ntechnology development to dissect genotype-phenotype associations and highlights how our\nframework provides an important foundation for such work.\nReviewer #2 (Remarks to the Author):\nIn this manuscript by Hsieh et al, the authors leverage a technique mitochondrial single-cell ATAC-\nseq, which they previously established to investigate heteroplasmy levels of specific variants from\npatients with mtDNA disease, (Nature Genetics 2022) to measure mtDNA mutational loads and\nsomatic mosaicism. To bench-mark their technique they use HEK-293 cell lines with a knock-in\nmutation in POLG which causes a proof-reading deficiency of the mtDNA polymerase and elevated\nmtDNA mutation rates. Interestingly in these lines, the variant allele frequencies (VAFs) of individual\ncells did not reach classical thresholds for metabolic dysfunction. However, when complex I and\ncomplex IV activity was measured in bulk cell cultures, there was metabolic impairment, suggesting\nmultiple mutations can act synergistically to affect mitochondrial function.\nThe authors then introduce two metrics: single-cell mtDNA mutations per million base pairs\n(scmtMPM) and heteroplasmy-weighted mitochondrial local constraint scores (scwMSS). These\nuseful measures allow comparative analysis of samples from subjects of different ages and disease\ngroups which is very useful for the research community. They then apply their benchmarked technique\nand measures to test datasets including control blood samples from young (5 year old) and older\n(47 year old) controls and 4 patients with the m.3243A>G mutation (aged 29-80). They show that the\nsomatic mtDNA mutation load per cell increases with age in control samples and suggest that the\nm.3243A>G mutation may affect the accumulation of somatic mtDNA mutations.\nSpecific comments:\n1. The authors have established a protocol to quantify mtDNA mutation load at the single cell level at\nan unprecedented resolution. Previous single cell work on somatic mtDNA mutations has been at the\nlevel of 100s of cells, this was 1000s. The data is very high resolution and supports previous studies\nwhich have shown similar mutational spectra and dynamics of age-related somatic mtDNA mutation\naccumulation.\nWe thank the Reviewer for the positive and encouraging feedback on the technical advances and\nresolution achieved by our methodology.\n2. This work is of high technical quality and provides a robust framework for assessing mtDNA\nload/constraints/evolution in other tissue/disease subgroups.\nWe thank the Reviewer for acknowledging the robustness and broad applicability of our framework\nfor assessing mtDNA mutational load and constraint.\n3. I think the title could be clearer, as although the technique is technically multi-omic, other than one\nexperiment towards the end where the nuclear ATACSeq data is used to subset immune-cells, this\npaper is really focussed on mtDNA sequence analysis. There is very little discussion of the impact (if\nthere is any) of cell type on somatic mtDNA mutagenesis. The majority of the analysis groups all of\nthe different cell types together. I would suggest either fully describing all of the data at the level of\ncell types, or just using the name of the technique in the title (mitochondrial single-cell ATAC-seq)\nwhich will show the focus of the paper.\nWe thank the Reviewer for this feedback regarding the extent to which we leverage the multi-omic\naspect of our dataset. We appreciate the Reviewer’s recognition that our primary focus is on high-\nthroughput single-cell mtDNA genotyping analysis, and the development of a broadly applicable\nframework for mtDNA mutational landscapes rather than an exhaustive characterization of nuclear\nscATAC-seq profiles. However, we note, that the scATAC-seq-derived nuclear data is the basis for\nall cell clustering and cell annotation analyses that is used throughout the manuscript. In addition, in\nthe revised manuscript, we have further expanded the multi-omic analysis specifically i) by comparing\nchromatin profiles of CTRL vs. POLGD274A cells (Extended Data Fig. 2), as well as ii) in the new\ngalactose treatment experiment, where we consider the integration of chromatin accessibility and\nmtDNA genotypes to be most informative (new Fig. 5 and Extended Data Fig. 8-9). Specifically, we\ndemonstrated mtDNA copy number amplification with concomitant chromatin accessibility changes in\ngenes related to mitochondrial stress under galactose-induced metabolic stress in the POLGD274A\nKI36 line. As such, we do consider it appropriate to retain ‘multi-omics’ in the title and hope for the\nReviewer’s agreement.\nNevertheless, the investigation of mtDNA mutational-burden–related chromatin accessibility changes\nin primary human cells represents an important and promising future direction, but we believe that a\nlarger cohort and substantial additional effort will be required to address this question rigorously,\nwhich we respectfully consider beyond the scope of the present study and its focus. However, moving\nforward, we believe our approach provides a framework for diverse future investigations to assess\nthe impact of somatic mtDNA mutagenesis across a variety of biological systems and contexts (e.g.,\naging and primary mitochondrial disorder), as we now also further outline in our revised Discussion.\n4. There appears to be three distinct populations of cells in the KI36 cell line, but the data is analysed\nas one population. Are there specific clonal variants that are driving the clustering in this line?\nWe thank the Reviewer for this comment. Yes, there are indeed various subclone-specific mtDNA\nvariants that are driving the sub-clustering of the KI36 line, select ones of which we visualized in\nFig. 4e. Please also consider the heatmap below, which shows 188 mtDNA variants that are shared\nand/or are specific to each subcluster (see Reviewer Fig. 4 below). Along these lines, we did note\nsubtle differences in mutational burden as evident by diverging mutational metrics for different gene\ngroups (Fig. 4c,d).\nReviewer Figure 4. Subclonal structure of KI36. (a) UMAP of the POLGD274A knock-in HEK293 clone KI36,\ncolored by clonotype ID (1-3), identified using the FindClonotype function in the R-package Signac. Manual\nsubcluster labels (C1–C3) are shown; each clonotype has a one-to-one match with a subcluster. (b) Heatmap\nshowing variant allele frequencies of 188 clonally informative mtDNA variants (rows) in individual KI36 cells\n(columns). Cells are grouped by clonotype (top color bar), highlighting the subclonal structure.\n5. Can you be more explicit when describing the confirmed pathogenic mutations on Mitomap? There\nare a number of different tables on there and I’m guessing you used the ‘confirmed’ table. It would be\ngood to state this to avoid confusion.\nWe thank the Reviewer for raising this important point for clarification. In the revised manuscript, we\nnow clarify that the table of confirmed pathogenic variants was obtained from the list of “Confirmed\nPathogenic Mutations” in MITOMAP. Variants were further filtered only for those annotated as\n“Confirmed (Cfrm)”, “Confirmed- Pathogenic (Cfrm [P])”, or “Confirmed- Likely Pathogenic (Cfrm [LP])”.\nThis is now described in the revised Methods section, and the list of all pathogenic variants evaluated\nin our study is now provided as a new Supplementary Table 4.\n6. For the scmtMPM score, from what I understand, this is normalised to coverage? Can you be\nexplicit that this is the total of the mtDNA aligned reads and not the overall sequencing depth of the\nATAC-Seq?\nWe thank the Reviewer for this question. Yes, the scmtMPM score is indeed normalized to the total\nnumber of bases that have been sequenced for the defined mtDNA region (mtDNA-aligned reads),\nrather than the total ATAC-seq read depth, which would also contain nuclear fragments. We have\nnow added this clarification to our revised Methods section to avoid any ambiguity.\n7. I don’t really understand how the conclusion was drawn that the 3243 mutation impedes somatic\nmtDNA mutation accumulation? In figures 5d and 5e both of these patients are relatively young (29\nand 35) and the total scmtMPM (-3243A>G) looks to be slightly lower than H47 (closest age-match).\nHowever, M60 and M80 are much higher supporting the accumulation of somatic mtDNA mutations\nwith age. What is the evidence that the 3243 variant impedes the accumulation of somatic variants?\nWhat happens if you draw similar plots for M60 and M80? I may be misunderstanding the figure, is\nthere an alternative way of presenting this data?\nWe apologize that this observation has not been sufficiently clearly described. We have deliberately\ndistinguished between the younger (M29 and M35) and older patients (M60 and M80), as at a more\nadvanced age, the pathogenic mt.3243A>G variant has been largely purified away. This is intimated\nin Fig. 6b as the cellular mutational burden significantly drops in the younger patients when excluding\nthe mt.3243A>G variant. In contrast, the burden remains largely unchanged in the older patients. The\ndecreasing levels of the mt.3243A>G variant with age are thereby consistent with prior\nobservations7,8,10,11. This may further suggest that somatic mtDNA mutation dynamics may indeed\nchange as a function of the declining mt.3243A>G variant over the lifetime of the patient. However,\nwe acknowledge the speculative nature of these interpretations.\nWith respect to the data shown in Fig. 6d,e, we have compared how the total scmtMPM changes\nwhen including the mt.3243A>G variant (panel d) vs. excluding it (panel e). What we noted is that in\ncells with a high mt.3243A>G variant allele frequency, the total scmtMPM score was almost entirely\ndriven by the presence of the pathogenic variant, given its removal leads to a significant decrease of\nthe captured mutational burden. As such, we interpreted that cells with a high pathogenic variant\nburden may be less tolerant to acquiring additional mutations. We also display the data as scatter\nplots in Extended Data Fig. 11e for the older patients (M60 and M80). However, the distribution does\nnot allow for a meaningful smooth kernel density representation, likely due to strong purifying selection\nand the resulting sparsity of cells at the high end of the mt.3243A>G VAF distribution. We understand\nthat there is a speculative nature to this interpretation, but we felt this is an intriguing observation\nrevealed through our single-cell somatic mtDNA burden analysis that we wanted to share with the\nreadership. It will be interesting to assess the extent of this phenomenon in additional patients,\nincluding those with mitochondriopathies due to other pathogenic mtDNA mutations. We have revised\nthe manuscript accordingly to better clarify the interpretations and future perspectives of the described\nobservation.\nReviewer #2 (Remarks on code availability):\nThis is not my area of expertise so I did not attempt to run the code.\nReviewer #4 (Remarks to the Author):\nThe study by Yu-Hsin Hsieh and colleagues presents a study in mitochondrial genomics through the\ndevelopment and application of novel single-cell metrics to quantify mitochondrial DNA (mtDNA)\nmutational burden and heteroplasmy dynamics. Utilizing mitochondrial single-cell ATAC-seq\n(mtscATAC-seq), the authors introduce two key metrics—single-cell mtDNA mutations per million\nbase pairs (scmtMPM) and heteroplasmy-weighted mitochondrial local constraint scores (scwMSS)—\nto capture genome-wide mtDNA variation at single-cell resolution. The study benchmarks these\nmetrics in POLGD274A knock-in cell lines and applies them to primary human immune cells from both\nhealthy donors and patients with mitochondrial encephalomyopathy (MELAS). While the work\naddresses an important challenge in mitochondrial genetics, offering potentially valuable insights into\nmitochondrial disease mechanisms, several critical issues and concerns warrant further attention.\nWe thank the Reviewer for their praise and the recognition of the importance of our work to the field\nof mitochondrial genetics.\nMajor Comments:\n1. The authors reported 455 variants in control cells versus 2484 and 4651 variants in the KI36 and\nKIA2 POLGD274A knock-in cells, respectively. Given the substantially higher sequencing depth in\nthe knock-in lines, the increased variant counts may largely reflect depth-related bias. Additionally,\nthe authors should also clarify how uneven mtDNA coverage across the mitochondrial genome was\naccounted for during variant calling and downstream analyses.\nWe thank the Reviewer for this important comment. In the revised manuscript, we have now\nreanalyzed 8x fold more deeply sequenced data from our HEK293 cell lines. This improved the raw\nreads from about 23,000 to 190,000 per cell (10x Genomics recommends up to 25,000 reads) and\nyielded a substantially increased per-cell mtDNA depth (median: CTRL = 55.6×, KI36 = 154.1×, KIA2\n= 255.2×). This further enabled us to perform downsampling analysis to systematically assess the\ninfluence of sequencing depth on mtDNA variant detection (see Reviewer Fig. 5 below; new\nExtended Data Fig. 7). As one may expect, the total number of detected variants increases with\nsequencing depth; however, the rate of increase began to plateau at approximately 50× mtDNA\ncoverage, even in the hypermutated POLGD274A knock-in background. Because mutation counts are\ninherently correlated with sequencing depth, we introduced the single-cell mtDNA mutations per\nmillion base pairs (scmtMPM) metric, which effectively normalizes for “position-wise” sequencing\ndepth by summing the total number of bases sequenced across the mtDNA regions of interest (i.e.,\nrespective gene, or all genes as part of a complex).\nRegarding the concern about uneven mtDNA coverage across the mitochondrial genome, our data\nshow near-uniform coverage (Fig. 1e), consistent with our previous demonstration and benchmarking\nof the mtscATAC-seq assay5. While some variation in coverage may be attributable to Tn5-related\ninsertion bias and other technical nuances, there are typically no significant differences in coverage\npatterns between samples.\nVariant calling is then performed using the standardized mgatk pipeline, which we have also\npreviously benchmarked in the same study by Lareau et al5. This pipeline incorporates a stringent\npaired-end read strand concordance threshold filter to ensure high-confidence variant detection.\nVariant allele frequencies (VAFs) are calculated based on position-specific coverage, and only sites\nsupported by at least five reads are retained for downstream analysis. Together, these measures\nminimize coverage-related artifacts and ensure reliable variant detection across the mitochondrial\ngenome.\nFurther, we note that differences in mtDNA coverage can also be biological, as different cell types\nhave different mtDNA copy numbers. In the context of the HEK293 lines, we note that the POLGD274A\nknock-in lines indeed have a higher copy number, which has previously been described as a\ncompensatory mechanism12–14.\nReviewer Fig. 5 (Manuscript Extended Data Fig. 7) – Impact of sequencing depth on mtDNA variant\ndetection in POLGD274Acell lines. (a) Linear regression of total mtDNA variant detection rate as a function of\nsequencing depths up to <40M reads. Data points represent the number of total mtDNA variants detected.\nRegression lines with 95% confidence intervals and equations are shown. (b) LOESS regression of total mtDNA\nvariant counts across the broader range of sequencing depths (10–200M reads). Data points represent the\nnumber of total mtDNA variants detected, and LOESS curves with 95% confidence intervals are shown. (c, d)\nViolin plots of mtDNA variant counts in single cells at (c) shallow sequencing depths (up to 40M reads) and (d)\nacross the broader range (up to 200M). The (c) linear regression and (d) LOESS modeling were performed on\nthe per-depth median number of detected variants. (e) Stream plot of relative heteroplasmy across the various\ndownsampled sequencing depths, depicting the overall variant composition. Each band in the stream plot thus\nrepresents the relative abundance, by normalizing to the sum summed pseudobulk VAFs of all variants at that\ndepth, of an individual mtDNA variant across downsampled sequencing depths. (f) Scatter plots show\npseudobulk VAFs of mtDNA variants or all detected variants in KIA2 when comparing datasets downsampled\nto 50M versus 200M total mtDNA reads. Each point represents a single mtDNA variant. Insets zoom in on\nvariants with VAF <1% to highlight concordance even at low allele frequencies. (g) Density plot of changes in\nestimated VAF (ΔVAF; 200M vs 50M) for all non-zero-VAF mtDNA variants (7,100 variants) across individual\nKIA2 cells (n = 1,894). (h) Violin plots of scmtMPM and scwMSS in single cells across the fixed subsampled\nmtDNA sequencing depth (10–200M mtDNA reads).\n2. The two metrics (scmtMPM and scwMSS) represent key contributions of this study. However,\nbiases in variant detection and variability in sequencing coverage could substantially influence their\naccuracy. The reported median mtDNA sequencing depth, ranging from approximately 10x to 40x,\nis very low for confident variant calling, particularly for low-frequency variants. Although the authors\nhave implemented measures to minimize these biases, further validation, such as targeted deep\nsequencing, is needed to confirm the accuracy of variant calls and to robustly support the conclusions\ndrawn from these metrics.\nWe thank the Reviewer for this comment, which also relates to comments #1 and #4 of Reviewer #1.\nAs suggested by the Reviewer and to assess the impact of sequencing depth on our mtDNA\nmutational burden metrics, we conducted 8x-fold deeper sequencing of our HEK293 dataset as also\nelaborated on in our response to the Reviewer’s prior comment. This enabled us to perform\ndownsampling analyses to simulate various sequencing depths across the different HEK293 cell lines.\nImportantly, both scmtMPM and scwMSS metrics remained stable across different downsampling\ndepths (new Extended Data Fig. 7h), demonstrating that our approach effectively mitigates depth-\nrelated bias in mutational burden quantification. Accordingly, we updated our HEK293 results with the\nmore deeply sequenced data in our revised manuscript (Figs. 1–4).\nWe further conducted downsampling analysis of the POLGD274A knock-in line data to assess the\ndetection of mtDNA variants as a function of their allelic frequencies at variable mtDNA sequencing\ndepth. As expected, this is very robust for higher allelic frequency variants, as shown in new Extended\nData Fig. 7. Notably, low-heteroplasmy variants (those <1% VAF) were also consistently recovered,\nand their pseudobulk VAF estimates remained highly concordant across depths (Extended Data\nFig. 7f). Moreover, single-cell mtDNA heteroplasmy estimates were highly consistent, with most per-\nvariant ΔVAF values within ±5% and a global ΔVAF distribution tightly centered at zero when\ncomparing up to 4-fold differences in sequencing depth (Extended Data Fig. 7g). Together, these\nanalyses suggest that allelic frequency estimates are largely stable once moderate coverage is\nachieved.\nFurther, we note that all HEK293 lines, including the control, would be expected to have the same\nbackground rate of potential false-positive mtDNA variants. As can be seen in Fig. 1d, the number of\nadditional variants detected in the POLGD274A knock-in line is, however, substantially larger than in\nthe control line. As such, it is extremely unlikely that this population and the large increase of detected\nvariants within is confounded by false-positive variant calls.\nRegarding the mtscATAC-seq and variants assay thereof, we emphasize that these have further been\nextensively benchmarked and applied2,4–8. Modified versions have further been adapted by members\nof the mitochondrial genetics community, such as for MitoPerturb-seq9. Earlier work further\ndemonstrated that single-cell heteroplasmies were well correlated across different single-cell genomic\nassays, including scRNA-seq and scATAC-seq4,5. In that work, we also demonstrated that pseudobulk\nheteroplasmies called from single-cell assays are highly correlated with bulk heteroplasmy\nmeasurements from bulk genomic sequencing. Please see, for example, Fig. 2F and Fig. S2E from\nLudwig et al., which we pasted below for the Reviewer’s consideration (Reviewer Fig. 2a,b).\n[editorial note: third party material redacted]\nWe also note our previously introduced paired-end, strand-aware variant calling workflow, as\nimplemented in the mgatk pipeline, provides higher specificity and more conservative filtering than\nstandard bulk mtDNA workflows5. As such, we would be uncertain to what extent bulk analysis,\nincluding via targeting of select variants among the pool of >10,000 detected variants, and challenges\naround sampling bias of very low frequency variants, can provide full orthogonal validation beyond\nour previously published evaluations2,4–8. However, we acknowledge the residual uncertainty\nsurrounding very low allele frequency heteroplasmic variants that may be only present in a small\nnumber of mtDNA copies or cells but would suggest that such variants are unlikely to meaningfully\nadd to a mutational burden analysis and/or be biologically impactful. We now discuss this in our\nrevised discussion, including the need for continued refinement of experimental and computational\nmethods, such as the integration of unique molecular identifiers and duplex-sequencing\napproaches15.\n3. The analysis assumes that variants unique to POLGD274A lines are de novo, but clonal expansions\nor bottlenecks could also influence variant distributions. Could the authors elaborate on how these\nfactors were considered or controlled for in their analyses?\nWe thank the Reviewer for highlighting this important aspect of mitochondrial genetics. We agree that\nstrictly “de novo” variants cannot be unequivocally established given the features of mtDNA\ninheritance (e.g., clonal expansion, bottlenecks) and technical factors (e.g., sampling bias,\nsequencing depth). Thus, we have revised our manuscript by replacing the term “de novo” with\n“POLGD274A-line specific” to reflect this concern accordingly. To mitigate potential confounding effects\nfrom clonal expansion and bottleneck effects, our downstream analyses only considered variants\nabsent from the parental HEK293 control cell line.\nWith respect to clonal expansions and bottleneck effects, these were not explicitly considered, as we\nwould be uncertain to what extent these may present potential confounders. Specifically, clonal\nexpansions would likely increase the VAF of the mtDNA variants and thus enhance their detectability.\nFor example, Fig. 4e shows a number of subclones specific to the line KI36. Given their recurrent\ndetection across a number of cells and their concentration in specific subclusters, one can consider\nthese very high-confidence variant calls.\nFurther, our read-downsampling analysis suggests the observed mutation events to be significantly\nmore prominent in the POLGD274A knock-in cell lines (Extended Data Fig. 7a-d), along with a lower\nVAF distribution (Fig. 2f,g). Together, this suggests that the variants are indeed likely to have\noccurred “de novo” in this genetic background. However, we agree that some residual uncertainty\nremains related to the confident detection of very low frequency mutations, and as such have\nreworded the description of these results and point out respective limitations in our revised\nmanuscript. Please see also our response to comment #2 regarding the background false-positive\nrate.\nWith respect to bottleneck effects, we would believe that longitudinal sampling experiments with a\nlarger population cohort may be more explicitly tailored to answer such questions, which will require\nmore comprehensive efforts, which we would respectfully consider to be more suitable for future\ninvestigations. Nonetheless, in our revised manuscript and as suggested by Reviewer #1, we have\nnow conducted a more defined pilot experiment in POLGD274A knock-in cell lines, where we utilized\ngalactose to create a metabolic bottleneck by enforcing OXPHOS metabolism (new Fig. 5). Notably,\nwe were able to detect various minor shifts in the cumulative heteroplasmy at the single cell level but\ndid not observe pronounced shifts in heteroplasmy for most variants in KI36 cells (Fig. 5e).\nInterestingly, a compensatory increase in mtDNA copy number was observed upon galactose\ntreatment, demonstrating that our approach is able to capture diverse nuances of mitochondrial\ngenetics, including small heteroplasmy shifts.\n4. The authors incorporate a mitochondrial genome constraint model based on population genetics\ndata from gnomAD. Have the authors considered the impact of population diversity, mitochondrial\nhaplotypes, as well as the cell-type-specific or context-dependent functional constraints?\nWe thank the Reviewer for raising this important question of human mitochondrial genetics. In our\nstudy, the mutational burden analyses were restricted to heteroplasmic variants. Haplogroup-\nassociated variants, which are typically present at homoplasmy, were excluded from the burden\nanalyses. In theory, as these variants are also typically invariant across the population if sampled\nfrom the same donor, these would also not necessarily present a confounding element to the burden\nanalysis.\nMoreover, in the original work describing the constraint model16, the authors did benchmark their MLC\nscore framework against large-scale population data (i.e., MITOMAP, ClinVar, HelixMTdb) to\ndemonstrate that associations between MLC-derived MSS and human phenotypes remained\nsignificant with or without adjusting for the haplogroup.\nWith respect to cell-type-specific or context-dependent functional constraints, we believe that no\nappropriate data would be available to definitely answer this question. The constraint model is\nultimately based on human population data that is most powered to assess primarily germline genetic\nvariation that is well captured at the level of human populations, but not at the level of individual cell\ntypes. Potentially, the Reviewer is alluding to the often cell-type-specific phenotypes observed in\npatients with mitochondriopathies. For example, our prior work demonstrated that pathogenic mtDNA\nvariants/deletions are not well tolerated in human subset of T cells7,8, which may be analogous to the\ncell-type specific effects the Reviewer may be referring to. However, to fully account for cell-type and\ncontext, we are likely to require the cell-/organ- and context-specific evaluation of the impact of\nvariants, which would present a substantial effort. While invaluable, we believe this is clearly beyond\nthe current scope of our manuscript.\n5. Have the authors explored integrating other variant annotations, like predicted pathogenicity scores\nor conservation, to help refine the MLC-based scoring?\nWe thank the Reviewer for this comment. Here, we would kindly refer the Reviewer to the original\narticle by Lake et al.16, where the MLC framework was systematically validated against population-\nlevel databases, including MITOMAP, ClinVar, and HelixMTdb. The study thereby demonstrated that\nMLC scores exhibit strong concordance with known pathogenicity and evolutionary conservation.\nAdditionally, the authors also provided experimental validation of the MLC score using base-editing\napproaches, further supporting its robustness. While we could envision further refinement of the MLC\nframework (e.g., inclusion of larger populations and of more diverse ancestry), we believe that our\ncurrent implementation already provides a strong genetic basis for our single-cell evaluations of\nmtDNA mutational burden. However, of course, the presented scwMSS framework is adaptable and\ncould readily integrate future refinements surrounding MLC scoring and potentially other future\nrelevant models, that may, for example, more comprehensively integrate data across diverse human\npopulations and more advanced variant annotation.\n6. Given that mtscATAC-seq simultaneously captures mtDNA variants and chromatin accessibility,\nhave the authors investigated potential interactions between mtDNA mutational burden and nuclear\nchromatin states? Such as, whether/how nuclear-mitochondrial interactions influencing variant\ndynamics? whether cells with higher mutational burden exhibit altered chromatin accessibility;\nAdditionally, single-cell Multiome technique, which profile both transcriptomes and chromatin\naccessibility in individual cells, could provide an even more powerful platform to investigate the\nfunctional consequences of mtDNA mutations. Have the authors considered applying or integrating\nsuch approaches?\nWe thank the Reviewer for this important comment. In our revised manuscript, we now more explicitly\nleverage the chromatin accessibility modality, to investigate nuclear gene regulatory responses to an\nelevated mtDNA mutational burden and galactose-induced metabolic bottlenecks. Specifically, we\nperformed a differentially accessible gene (DAG) analysis comparing i)POLGD274A knock-in vs. CTRL\nHEK293 cells (new Extended Data Fig. 2), and ii) POLGD274A knock-in KI36 line under glucose vs.\ngalactose conditions (new Fig. 5). The first analysis revealed profound chromatin remodeling in both\nKI36 and KIA2 lines. Shared upregulated and downregulated DAGs were enriched for genes involved\nin extracellular matrix remodeling (e.g., protocadherins gene family, PCDHs), transcriptional\nregulation (e.g., zinc-finger transcription factors, ZNFs), mtDNA maintenance (e.g., MGME1), and\nrespiratory chain assembly, indicative of broad nuclear epigenetic responses to an elevated mtDNA\nmutational burden. Similarly, DAG analysis comparing glucose- vs. galactose-treated POLGD274A cells\n(KI36 line) identified genes involved in mitochondrial/ER stress and inflammatory signaling, indicative\nof epigenetic rewiring and adaptation in the presence of a profound metabolic bottleneck. Together,\nthese findings showcase the capability of mtscATAC-seq to capture biologically meaningful nuclear-\nmitochondrial interactions linked to mitochondrial dysfunction.\nWith respect to the Reviewer’s suggestion regarding the utility of the single-cell Multiome approach,\nwe agree that this approach is a powerful platform to investigate the functional consequences of\nmtDNA mutations and is highly complementary to our current mtscATAC-seq framework. Along these\nlines, our prior work identified extensive gene regulatory changes in human primary immune cells with\npathogenic mtDNA deletions revealed via scRNA-seq8. However, the number of samples were still\nlimited. To meaningfully extend this line of work, we are aiming to profile larger cohorts of healthy\ndonors as well as mitochondriopathy patients to be sufficiently well-powered and account for nuclear\nas well as mitochondrial genetic heterogeneity. Here, and for the Reviewer’s consideration, we have\nincluded some of our preliminary work (see Reviewer Fig. 6 below), demonstrating the potential for\nfuture applications of our MMB metrics.\nSpecifically, we stratified cells with a very low vs high mutational burden in a select healthy donor and\na MELAS patient (Reviewer Fig. 6a,c), and performed differential transcriptome analysis in a cell\ntype-specific manner. While a substantial number of differentially expressed genes were identified in\nboth the healthy donor and the MELAS patient cells in various cell types, the results yet lacked\nstatistical significance (Reviewer Fig. 6b,d). Additionally, in the MELAS context, the differences could\npotentially be dominated by the strongly pathogenic mt.3243A>G variant. While the investigation of\nmtDNA mutational-burden–related chromatin accessibility and transcriptome changes in primary\nhuman blood cells represents an important and promising future direction, we believe that for these\ntypes of analysis, one would require profiling of a significantly larger number of cells and donors,\nanalogous to eQTL (expression quantitative trait locus) studies, which aim to identify how genetic\nvariations influence gene expression levels by mapping genetic variants to gene expression\nphenotypes17–19. However, we hope the Reviewer agrees with our reasoning, that these would readily\npresent independent and stand-alone efforts. Importantly, and moving forward, we believe our\napproach provides a framework for a variety of future investigations to assess the impact of somatic\nmtDNA mutagenesis across a variety of biological systems and contexts (e.g., aging and primary\nmitochondrial disorders).\nReviewer Fig. 6 Preliminary integrated analysis of transcriptomic profiles and mtDNA mutational burden\nin a healthy donor and a MELAS patient PBMC datasets. (a) Joint transcriptome–chromatin accessibility co-\nembedded weighted-nearest-neighbor (WNN) UMAP projection of PBMCs from a healthy donor, colored by\nmajor immune cell types (left) and scmtMPM (right). (b) Violin plots of selected differentially expressed\ntranscripts comparing high-scmtMPM cells (red; scmtMPM > 65th percentile) and low-scmtMPM cells (blue;\nscmtMPM < 35th percentile), stratified by cell type. (c) WNN-UMAP projection of PBMCs from a MELAS patient,\ncolored by major immune cell types (left) and scmtMPM (middle), and the VAF of the pathogenic variant\nmt.3243A>G. (d) Violin plots of selected differentially expressed transcripts in the MELAS patient, using the\nsame high- vs low-scmtMPM comparison as in (b). Differential expression analysis was performed using the\n`FindMarkers` function in R-package Seurat (v4.0.3) using the default Wilcoxon Rank Sum test.\n7. It would strengthen the impact of the study if the authors included functional assays to see how\nmtDNA variants actually affect cells. Further linking the scmtMPM and scwMSS metrics to cellular\nphenotypes would add important biological context to the current study.\nWe thank the Reviewer for this important suggestion. As also suggested by Reviewer 1 (comment\n#5), we have now included new data where the POLGD274A knock-in cell lines were cultured in the\npresence of galactose, enforcing OXPHOS utilization. This enabled us to evaluate potential selection\ndynamics of cells and mutations under a metabolic bottleneck, thereby providing an additional layer\nof phenotypic characterization. We present these results in new Fig. 5 and Extended Data Figs. 8-\n9. Specifically, the POLGD274A knock-in cell lines exhibited a strong proliferation defect in galactose-\ncontaining medium, indicating a pronounced mitochondrial defect. Notably, instead of pronounced\nshifts in heteroplasmy, we observed mtDNA amplification with concomitant chromatin accessibility\nchanges in genes related to mitochondrial stress, suggesting the primary modes of compensation in\nthis experimental setting.\nHowever, as we now note in the revised Discussion, we acknowledge that varying the concentration\nor duration of treatment could result in different mutation dynamics. Moreover, the relatively uniform\nmutational burden observed across the KI36 line suggests that galactose treatment may have affected\nall cells similarly, rather than selecting for a more tolerant (clonal) subpopulation. While we consider\nfurther studies in this direction to be promising, we respectfully consider them as best suited for future\ninvestigations outside of the scope of this manuscript. Nevertheless, the decreased complex I and IV\nactivities (Extended Data Fig. 1a) and substantial galactose-induced cell death/growth defects\nobserved in the POLGD274A knock-in cell lines (Extended Data Fig. 8a) provides strong additional\nevidence that the elevated mutational burden significantly compromises mitochondrial function.\nMoving forward, and as we now mention in our revised Discussion, the direct integration of\nalternative phenotyping approaches, such as single-cell metabolite profiling20 in conjunction with\nscmtMPM/scwMSS measurements are likely to provide a more direct evaluation on how mtDNA\nmutational burden affects cells. However, we are currently unaware of an approach that would enable\nready integration of the detection of both modalities. Future efforts will be required to catalyze such\nmulti-omic analyses beyond the level of classical genomic approaches.\n8. Could the authors expand on the translational implications of their findings. In particular, discussing\nhow feasible it might be to apply the scmtMPM and scwMSS metrics in clinical or therapeutic settings,\nalong with any potential challenges, would add valuable perspective.\nWe thank the Reviewer for this valuable suggestion. In the revised manuscript, we have now added\na new paragraph in the Discussion outlining the potential translational applications of our scmtMMB\nmetrics. Specifically, we highlight recent discoveries linking mtDNA mutations and heteroplasmy to\nclinical outcomes in oncology and hematological malignancies. Moreover, elevations in mutational\nlevels have been associated with aging phenotypes in mutator mice21–23. As such, our single-cell\nframework is expected to be broadly applicable across different biological contexts, though\nappropriate cohort sizes will be necessary to be adequately powered to detect respective\nassociations. Along these lines, base-editing approaches are likely to catalyze the experimental\nvalidation of such observations.\n9. The manuscript currently lacks enough detail about the statistical analyses used. The authors\nshould include clear descriptions of the statistical methods applied throughout the study to better\nsupport their findings and improve transparency.\nWe thank the Reviewer for this important comment. We have accordingly revisited our descriptions\nof the applied methodology and statistical approaches to add more detail and clarity.\nOther Comments:\nIt would be helpful if the authors provided a summary table listing all variants detected in this study,\nincluding average sequencing depth, the number of cells carrying each mutation within each group,\nand other relevant information.\nWe thank the Reviewer for this suggestion. We now provide a comprehensive Supplementary\nTable 1 listing all high-confidence variants, coverage, cell counts, heteroplasmy, and annotations.\nReferences\n1. Stoeckius, M. et al. Simultaneous epitope and transcriptome measurement in single cells. Nat\nMethods 14, 865–868 (2017).\n2. Mimitou, E. P. et al. Scalable, multimodal profiling of chromatin accessibility, gene expression\nand protein levels in single cells. Nat Biotechnol 39, 1246–1258 (2021).\n3. Peeva, V. et al. Linear mitochondrial DNA is rapidly degraded by components of the replication\nmachinery. Nat Commun9, 1727 (2018).\n4. Ludwig, L. S. et al. Lineage Tracing in Humans Enabled by Mitochondrial Mutations and Single-\nCell Genomics. Cell 176, 1325–1339.e22 (2019).\n5. Lareau, C. A. et al. Massively parallel single-cell mitochondrial DNA genotyping and chromatin\nprofiling. Nat Biotechnol39, 451–461 (2021).\n6. Lareau, C. A. et al. Mitochondrial single-cell ATAC-seq for high-throughput multi-omic detection\nof mitochondrial genotypes and chromatin accessibility. Nat Protoc18, 1416–1440 (2023).\n7. Walker, M. A. et al. Purifying Selection against Pathogenic Mitochondrial DNA in Human T\nCells. N Engl J Med383, 1556–1563 (2020).\n8. Lareau, C. A. et al. Single-cell multi-omics of mitochondrial DNA disorders reveals dynamics of\npurifying selection across human immune cells. Nat Genet 55, 1198–1209 (2023).\n9. Burr, S. P. et al. MitoPerturb-Seq identifies common and gene-specific single-cell responses to\nmitochondrial DNA depletion and heteroplasmy. bioRxiv doi:10.1101/2025.07.08.663208.\n10. Franklin, I. G. et al. T cell differentiation drives the negative selection of pathogenic\nmitochondrial DNA variants. Life Sci Alliance 6, (2023).\n11. Zhang, J. et al. Antigen receptor stimulation induces purifying selection against pathogenic\nmitochondrial tRNA mutations. JCI Insight 8, (2023).\n12. Hedberg-Oldfors, C. et al. Deep sequencing of mitochondrial DNA and characterization of a\nnovel mutation in a patient with arPEO. Neurol Genet 6, e391 (2020).\n13. Salminen, T. S. et al. Mitochondrial genotype modulates mtDNA copy number and organismal\nphenotype in Drosophila. Mitochondrion 34, 75–83 (2017).\n14. Filograna, R. et al. Modulation of mtDNA copy number ameliorates the pathological\nconsequences of a heteroplasmic mtDNA mutation in the mouse. Sci Adv 5, eaav9824 (2019).\n15. Lawson, A. R. J. et al. Somatic mutation and selection at population scale. Nature647, 411–\n420 (2025).\n16. Lake, N. J. et al. Quantifying constraint in the human mitochondrial genome. Nature 635, 390–\n397 (2024).\n17. Nathan, A. et al. Single-cell eQTL models reveal dynamic T cell state dependence of disease\nloci. Nature606, 120–128 (2022).\n18. Zhang, Z. et al. Unveiling genetic signatures of immune response in immune-related diseases\nthrough single-cell eQTL analysis across diverse conditions. Nature Communications 16, 7134\n(2025).\n19. Hong, S. E. et al. Single-cell eQTL analysis identifies genetic variation underlying metabolic\ndysfunction-associated steatohepatitis. Nature Genetics 57, 1638–1648 (2025).\n20. Delafiori, J. et al. HT SpaceM: A high-throughput and reproducible method for small-molecule\nsingle-cell metabolomics. Cell 188, 6028–6043.e11 (2025).\n21. Trifunovic, A. et al. Premature ageing in mice expressing defective mitochondrial DNA\npolymerase. Nature429, 417–423 (2004).\n22. Kujoth, G. C. et al. Mitochondrial DNA mutations, oxidative stress, and apoptosis in mammalian\naging. Science 309, 481–484 (2005).\n23. Maclaine, K. D., Stebbings, K. A., Llano, D. A. & Havird, J. C. The mtDNA mutation spectrum in\nthe PolG mutator mouse reveals germline and somatic selection. BMC Genom Data 22, 52\n(2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for their thorough and meticulous response to my first review, and particularly for the generation of new\ndata. This, of course, raises further questions which they may want to consider.\n1. The authors generated more reads from the three cell lines. Does the increased read depth influence their results and\nconclusions?\n2. Why is the read depth different between the cell lines? I presume this reflects a difference in mtDNA content. Has this\nbeen validated using an orthogonal method?\n3. The comparison of heteroplasmy value using pseudo and real-bulk analysis really mises the point I was making – that low\ncoverage compromises the ability to call and measure heteroplasmy reliably. I do not think the current data really addresses\nthis.\n4. The authors should reword their section on the ‘metabolic bottleneck’. As written, it conflates two issues: (1) a bottleneck –\nif present you should see heteroplasmy levels changing by accelerated random drift – or the random loss/gain of\nheteroplasmies; (2) the metabolic stress which they suggest has killed a large proportion of cells. If this is the case, then you\nwould expect less pathogenic variants after the metabolic challenge. It is not clear to me whether the authors saw this or not.", "reviewer_2": "I enjoyed reading the revised version of this manuscript. I am happy that that authors have addressed my comments and the\nadditional analyses have strengthened the manuscript.", "reviewer_3": "I thank the authors for their thorough responses. My previous concerns have been addressed, and I have no further\ncomments on the manuscript."}, "author_response": "We thank the Reviewers for their positive feedback, which further helped to improve the manuscript. We provide a\ndetailed point-by-point response to all comments in blue and have modified the manuscript based on the valuable\nsuggestions and journal formatting guidelines, including shortening of the main text.\nDetailed Response to Reviewers\nReviewer #1 (Remarks to the Author):\nI thank the authors for their thorough and meticulous response to my first review, and particularly for the generation of\nnew data. This, of course, raises further questions which they may want to consider.\n1. The authors generated more reads from the three cell lines. Does the increased read depth influence their results\nand conclusions?\nWe thank the Reviewer for this follow-up question. The increased sequencing depth improved the sensitivity and\nrobustness of mtDNA variant detection, particularly for low-frequency variants. However, it did not change the overall\nresults or conclusions of our study. While deeper sequencing detected additional variants, the relative differences\nbetween CTRL and POLGD274A knock-in lines remained comparable, and the distributions of scmtMPM and scwMSS\nwere stable across a wide range of sequencing depths, as illustrated in our downsampling analyses (Supplementary\nFig. 7).\n2. Why is the read depth different between the cell lines? I presume this reflects a difference in mtDNA content. Has\nbeen validated using an orthogonal method?\nWe agree that the observed differences in read depth likely reflect biological differences in mtDNA content and the\naccumulation of linearized mtDNA, as previously described using an orthogonal method such as Southern blot in\nPeeva et al. (2018, PMID: 29712893) and Trombly et al. (2023, PMID: 37237953). More broadly, our prior work has\nshown that the percentage of recovered mitochondrial reads varies across tissues, as revealed by our reanalysis of\nGenotype-Tissue Expression (GTEX) data. For example, this analysis showed particularly high mitochondrial read\npercentages and coverage in human muscle, brain, and heart compared with other tissues, which aligns well with\ntheir higher mitochondrial content. Please see Figure 4a-c from DOI: 10.1016/j.cell.2019.01.022), which we have also\ncopied here for the Reviewer’s consideration.\n[editorial note: third party material redacted]\n3. The comparison of heteroplasmy value using pseudo and real-bulk analysis really mises the point I was making –\nthat low coverage compromises the ability to call and measure heteroplasmy reliably. I do not think the current data\nreally addresses this.\nWe thank the Reviewer for raising this concern again. The reliable detection of very low-frequency variants is\ngenerally challenging in genomics, whether in nuclear or mitochondrial genome sequencing. We believe this may\nultimately be addressable only via very high-fidelity sequencing platforms, such as NanoSeq, a duplex sequencing\nmethod with an error rate lower than five errors per billion base pairs, as introduced by Lawson et al. (2025, PMID:\n41062696). However, these technical principles have yet to be adopted for single-cell sequencing. However, as\nexplained in the prior round of review, the use of strand-concordance already provides a quite conservative threshold\nfor variant calling in our single-cell data. Nevertheless, we acknowledge in the revised discussion that the detection of\nvery low-frequency variants must be interpreted with caution, both technically and biologically, with respect to their\nfunctional impact.\n4. The authors should reword their section on the ‘metabolic bottleneck’. As written, it conflates two issues: (1) a\nbottleneck – if present you should see heteroplasmy levels changing by accelerated random drift – or the random\nloss/gain of heteroplasmies; (2) the metabolic stress which they suggest has killed a large proportion of cells. If this is\nthe case, then you would expect less pathogenic variants after the metabolic challenge. It is not clear to me whether\nthe authors saw this or not.\nWe thank the Reviewer for this important conceptual clarification and agree that the original wording conflated two\ndistinct processes. We have accordingly revised the manuscript for clarity. As demonstrated in our original\nsubmission, pathogenic variants were less prevalent than other variant types, which we interpreted as cells with\npathogenic variants being selected against over the course of culture, independent of galactose treatment. As such,\n1\nwe believe a different experimental setup is likely necessary to reaffirm the Reviewer’s notion, as presented in the\ndiscussion of our work.\nReviewer #2 (Remarks to the Author):\nI enjoyed reading the revised version of this manuscript. I am happy that that authors have addressed my comments\nand the additional analyses have strengthened the manuscript.\nWe sincerely thank the Reviewer for this praise and the valuable feedback during the review of our work.\nReviewer #3 (Remarks to the Author):\nI thank the authors for their thorough responses. My previous concerns have been addressed, and I have no further\ncomments on the manuscript.\nWe are grateful to the Reviewer for their time and suggestions to improve our manuscript.\n2"}]

Decision letter