Single-cell lineage tracing by endogenous mutations enriched in transposase accessible mitochondrial DNA

1Resolving lineage relationships between cells is necessary to understand the fundamental mechanisms underlying normal development and the progression of disease. In recent years, new methods have emerged to enable cell lineage tracking with increasing resolution, leading to substantial biological insights. Specifically, genome editing of reporter constructs via CRISPR-Cas9 allowed synthetic reconstruction of cell lineage relationships in model organisms, and has been coupled with transcriptome profiling to inform cell fates. These prospective ‘mutate-and-record’ methods provide powerful tools to resolve the developmental origin of cells in genetically engineered cells and organisms, but cannot be utilized in living humans, archival clinical samples, or any wild type organism. Given these limitations, retrospective lineage tracing using endogenous genetic markers is an alternative solution. Recent advances in sequencing enable naturally occurring somatic mutations to be used as lineage markers, which usually required single-cell genome sequencing to capture the sparse genetic information. Regions with high mutation rates, such as microsatellite repeats, retrotransposons, and copy-number variants, has been used to resolve the lineage relationship for normal or cancerous tissue samples. These methods reduce the cost of whole genome sequencing, but still lack information on cell phenotypes.

2Simultaneous measurement of the lineage relationship and cell fates is ultimately required to address many biomedical questions. Here we describe EMBLEM (Epigenome and Mitochondrial Barcode of Lineage from Endogenous Mutations), a strategy to track cell lineage using endogenous mitochondrial DNA variants in ATAC-seq data. The end result of EMBLEM is single-cell lineage information and rich global epigenomic profile from the same individual cells (Figure 1A and Figure 1—figure supplement 1).

3We illustrate the utility of EMBLEM in human blood progenitor cells to clarify the process of pre-leukemic clonal evolution and the emerging biology of clonal hematopoiesis.

4Assay of Transposase-Accessible Chromatin by sequencing (ATAC-seq) is a sensitive method used to study chromatin accessibility profiles in diverse cell types and organisms. During DNA transposition and amplification in cells, mitochondrial DNA is also amplified at the same time (Figure 1A). Mitochondrial DNA (mtDNA) is a ~ 16 kb circular genome with ~10 fold higher mutation rate compared to the nuclear genome. Hence, mtDNA incrementally accumulates unique, irreversible genetic mutations that are passed on to daughter cells even in healthy humans and may be used for lineage tracing,. The majority of somatic mtDNA mutations are noncoding and thought to be passenger. Importantly, the number of mitochondria (and therefore mtDNA) range from several hundreds to >10,000 per cell in different cell types, facilitating robust mtDNA analysis even from a single cell.

5We first observed that ATAC-seq effectively enriches for mtDNA. While mtDNA is present in many kinds of DNA sequence libraries, it is substantially enriched in ATAC-seq libraries due to the fact that mtDNA is not chromatinized and is therefore highly accessible (Supplementary file 1). ATAC-seq enables a 17-fold or greater enrichment of mtDNA compared to exome sequencing or whole genome sequencing in GM12878 human B cells (Figure 1B), leading to an average ~18,000X coverage of mtDNA (Figure 1—figure supplement 2A). With this coverage, we detected 27 mitochondrial variants from GM12878 cells (Figure 1C). 13 of these variants have a variant allele frequency (VAF) greater than 90%, which are known as homoplasmic variants (Figure 1—figure supplement 2B). We also detected 14 low frequency mitochondrial DNA variants, with VAFs ranging from 0.1% to 24% (Figure 1C and Figure 1—figure supplement 2C). Similar results for mtDNA enrichment were observed in human K562 cells (Figure 1—figure supplement 3, Supplementary file 1).

6The VAF from bulk ATAC-seq data represents the average of the allele frequencies of the cell population. A 25% VAF may arise from 25% of cells in the population with a homoplasmic variant, or alternatively arise from 100% of cells all having a quarter of their mitochondria with the variant allele (Figure 1D). To distinguish between these two models, we analyzed single-cell ATAC-seq data from GM12878. For four mtDNA variants (VAF between 0.5% ~ 24% at population level), we find that a mixture of both models is in action for different variants (Figure 1E). For instance, mtDNA mutation 3082 is widely spread among single cells, but at low frequency per cell. Because it is extremely unlikely (see Materials and method) that the identical mutation arose independently in every single cell, cells sharing the same mitochondrial mutations are inferred to have descended from the same ancestral cell. These results suggest that even low frequency heteroplasmic mtDNA mutations can be exploited for lineage tracing.

7To prove the principle that somatic mitochondrial mutations can track cells from the same ancestor and to quantify the power of lineage mapping, we next applied EMBLEM to primary blood cells from patients with acute myeloid leukemia (AML). Human AML is organized as a hierarchy: a hematopoietic stem cell first acquires an initiating mutation in one of a number of chromatin modifier genes, previously termed as ‘pre-leukemic’ hematopoietic stem cell (pHSC). pHSCs are functionally normal and are not able to transplant AML, but upon accumulation of additional mutations, they give rise to leukemic stem cells (LSCs) that are able to self-renew and recapitulate AML disease upon transplantation. Finally, LSCs give rise to the bulk leukemic blast cells in AML. Targeted exome sequencing in these samples have identified somatic mutations in tumor suppressor genes and oncogenes that link the lineage relationship of pHSCs, LSCs and blasts, providing the ground truth for our analyses.

8We applied EMBLEM to the ATAC-seq profiles of FACS-purified LSCs and leukemic blasts first (Supplementary file 1). Using high-confidence mtDNA mutations, detected both from bulk ATAC-seq and single cell ATAC-seq, we found the LSC and blast populations not only shared the same heteroplasmic variants, but also showed similar distribution and allele frequency at the cellular level (Figure 1—figure supplement 4). These results indicate the two populations are identical at the genetic level, but divergent at the epigenomic level, consistent with previous studies. In patient SU353, we identified four diagnostic mtDNA mutations in the same cell (Figure 1F), which indicates these four mitochondrial variants already co-existed in the ancestral cell (see Materials and method). With the assumption that all these LSCs and blasts are clonal, we further quantified the detection rate of each mtDNA variant as a function of allele frequency and sequencing depth (Figure 1G). We found that when a single variant allele has a frequency greater than 20%, the detection rate can be up to 90% with >20X coverage (e.g. site 6776). In contrast, when the variant allele has a frequency lower than 1%, the detection rate drops to 20% when the coverage is below 100X (e.g. site 6705). While high drop-out rate is a common challenge for single-cell technologies, computational imputation of the missing information from single cell data can address this problem. When multiple mtDNA variants are co-detected in multiple single cells, we can infer their origin and linkage in the ancestral cell (see Materials and method). Thus, cells containing any one of these variants will still inform their origin from the same lineage. With any combination of the four variants, 90% (sensitivity) of the cells can be unambiguously assigned to the correct lineage with just 20x mtDNA coverage (Figure 1H). Furthermore, two mtDNA mutations identified in other cells (e.g. pHSC specific site 2967,6268) were never detected (false positive = 0) in LSCs and blasts (Figure 1—figure supplement 4), showing a high specificity of the method. Similar performance of single cell lineage tracing for another patient (SU070) is shown in Figure 1—figure supplement 5. These results demonstrate that somatic DNA mutations in the mitochondrial genome are a powerful endogenous marker to identify clonal cell populations.

9To expand on these findings to additional different cell lineages, we applied EMBLEM to bulk ATAC-seq data from sorted blood cells from healthy human donors and patients with AML(Supplementary file 1). We identified heteroplasmic mtDNA mutations in multiple cell populations of primary blood cells from healthy donors and all AML patients (Figure 2—figure supplement 1A, Supplementary file 2). The heteroplasmic mtDNA mutations showed a similar mutant spectrum as observed by previous studies using cancer genomic data (Figure 2—figure supplement 1B and C).

10Furthermore, EMBLEM, not only confirmed the previous lineage hierarchy of AML, but also extended the previous model of pHSC heterogeneity (Figure 2A). In the AML cases with LSCs sequenced by ATAC-seq, the LSCs and their corresponding leukemic blasts have nearly identical heteroplasmic mtDNA mutations (Figure 2B–C and Figure 2—figure supplement 1D), suggesting a direct lineage relationship and short generation history between LSCs and blasts. We then examined whether any of the mtDNA variants present in LSCs can be seen in the pHSCs, where the first leukemia-associated protein-coding mutations have already occurred in functional normal hematopoietic stem cells. We detected blast-associated mtDNA mutations in pHSCs in all 11 cases. Interestingly, we also detected additional heteroplasmic mtDNA mutations present specifically in pHSCs (Figure 2C). In the 11 cases we investigated, seven cases have pHSC-unique heteroplasmic mtDNA mutations (Figure 2—figure supplement 1D and E), a previously unrecognized level of pHSC heterogeneity. pHSCs are capable of long-term self-renewal and possess a clonal growth advantage, allowing them to clonally outcompete normal HSCs. Indeed, the clonal frequency of pHSCs is a poor prognostic factor for overall survival in AML. Our discovery of pHSCs with distinct heteroplasmic mtDNA mutations suggests the existence of multiple distinct sub-clones of pHSCs in AML patients.

11To validate the heterogeneity of pHSCs inferred from EMBLEM of bulk cell populations, we performed single-cell ATAC-seq of HSCs from AML patient SU353, which exhibited both a high burden of pre-leukemic somatic coding gene mutations and high frequency of pHSC-specific heteroplasmic mtDNA mutations. We identified the heteroplasmic mtDNA variants from each single cell, which separated the HSCs into three lineages: Two clonal subpopulations termed ‘clone 1’ (18 cells) and ‘clone 2’ (104 cells), and a third population with no mtDNA variants despite sufficient mtDNA coverage (pHSC with WT mtDNA, 31 cells) (Figure 2D). Notably, clone 2 possessed pHSC-specific mtDNA mutations, while clone 1 possessed mtDNA mutations shared with LSCs, indicating clone1 is the lineage precursor of AML. These results confirm that multiple pHSC clones arise in AML patients, and one subclone eventually evolved to become the LSC (Figure 2—figure supplement 2A).

12Finally, we related the clonotype of pHSCs to their single-cell chromatin accessibility profiles. We interrogated the patterns of active DNA elements and enriched transcription factor motifs in sequential stages of AML development from the same patient, and contrasted with HSCs from normal donors using ChromVAR (Figure 2E and Figure 2—figure supplement 2B). The chromatin accessibility profiles of pHSCs are more similar to HSCs than to LSCs or leukemic blasts. The greatest deviation between HSC and other cell types occurred at DNA binding motifs of the transcription factor Jun/Fos, a known key regulator of HSC biology (Figure 2F). Furthermore, the three lineages of pHSCs revealed by mtDNA mutations also showed distinctive chromatin profiles (Figure 2G). Clone 1 pHSC, which gives rise to the LSC and AML leukemia, is already more similar to LSCs and blasts in its chromatin accessibility. In contrast, clone 2 that comprises the larger fraction of pHSCs exhibited variable chromatin profiles at the single-cell level that spanned the range of normal HSCs, pHSC with WT mtDNA(WT cells), is also diverged from normal HSCs. Thus, both lineage tracing and single cell epigenomic states indicate clone 1 as the original stem cell of the AML in patient SU353. Supervised comparison of the chromatin accessibility profiles among these clonal sub-populations further identified distinct and significantly enriched transcription factor motifs (Figure 2H and Figure 2—figure supplement 2C–E). These results indicate the heterogeneity of HSCs from AML patients both on a genetic and epigenomic level.

13We present a computational strategy to combine cell lineages tracing by endogenous mtDNA mutations and chromatin accessibility profiling in the same cell using single-cell ATAC-seq data. This approach is applicable to any eukaryote, does not require genetic engineering or genome editing, and is cost effective as the lineage information comes ‘for free’ on top of epigenomic insights. The relative merits of mtDNA vs other genetic markers for lineage tracing are outlined in Supplementary file 3. An important advantage of EMBLEM is that we enable clonotype tracing in existing ATAC-seq data sets and hierarchical lineage construction from ATAC-seq that thousands of labs have already generated. All future ATAC-seq data acquired for other inquiries will also have the benefit of lineage information. EMBLEM may also be extended to other single cell technologies, in which mtDNA is sequenced. We show that EMBLEM is successful even with low frequency heteroplasmic mutations, detection of rare clones in a population, and authentic clinical samples. With advances in the throughput and depth of single-cell genomic technologies, we believe EMBLEM may be a powerful tool to bring insight for many biomedical questions, including development, regeneration, immunity, and cancer with integration of genotype and phenotype information from the same cell. During revision of this work, reported the feasibility of using mtDNA and single cell genomics for lineage tracing, which independently validates the potentially broad utility of this approach.

14Although powerful and broadly applicable, mtDNA lineage tracing also has its limitations. One limitation of this method is absence of mtDNA mutations in cells and tissues of embryos and young animals, which precluded us from applying EMBLEM to published scATAC-seq data of early animal development. Moreover, the possibility of selective mitochondrial inheritance or intercellular mitochondria transfer may affect the accuracy of inferred lineages. On the other hand, asymmetric transmission of mitochondria would not necessarily affect cellular lineage tracing, as long as the variant alleles are randomly segregated. Using scATAC-seq data from a mixing experiment with human and mouse cells, we found species-specific mtDNA always paired with species-specific nuclear genomic DNA (Figure 2—figure supplement 4). These results suggest that mitochondrial horizontal transfer is not a confounder of our study and does not universally occur between cells. The aforementioned two scenarios reflect the potential uncoupling of nuclear and mitochondrial genomes, which would be of interest to investigate by EMBLEM in combination with other nucDNA tracing methods.

15mtDNA lineage tracing produced new insights concerning the pHSC, the human hematopoietic stem cell that suffers the first oncogenic mutation in AML evolution. Our results add to the evidence that the pHSC population is heterogeneous, with evidence of multiple mtDNA clones. Unexpectedly, the pHSC lineage that gives rise to the subsequent acute myeloid leukemia is not the lineage with the best competitive potential among pHSCs, as the leukemogenic lineage is often in the minority. pHSC burden is a strong poor prognostic predictor of AML survival. It is widely believed that the association between high pHSC burden and poor AML patient prognosis reflected the enhanced self-renewal and competitive ability of the mutant pHSC. Our analysis suggests that high pHSC burden may reflect the diversity of pHSCs or the underlying mutational processes. These alternative interpretations of the link between pHSC burden and poor clinical prognosis should be addressed in future studies.

16Aligned bam files for GM12878 whole exome, low coverage whole genome, and PCR free whole genome sequence, were downloaded through phase 3 release of 1000 genomes (ftp://ftp.1000genomes.ebi.ac.uk).

17The alignment files were accessed via the following ftp links:

18ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA12878/exome_alignment/NA12878.mapped.ILLUMINA.bwa.CEU.exome.20121211.bam

19ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA12878/alignment/NA12878.mapped.ILLUMINA.bwa.CEU.low_coverage.20121211.bam

20ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA12878/high_coverage_alignment/NA12878.mapped.ILLUMINA.bwa.CEU.high_coverage_pcr_free.20130906.bam

21ATAC-seq and single cell ATAC-seq data for GM12878 generated by Buenrostro et al. were downloaded through GEO with accession number GSE47753 and GSE65360, respectively. Bulk ATAC-seq data from normal donors and AML patients generated by, were downloaded through GEO with accession number GSE74912. Single cell ATAC-seq data for leukemia stem cell and leukemic blasts generated by the same study were downloaded through GEO with accession number GSE74310. Single cell ATAC-seq from normal HSC generated by, were downloaded through GEO with accession number GSE96772. Supplementary file 1 summarized the detail information of all the datasets used in this study.

22Sequencing reads from ATAC-seq were aligned to the reference genome by BWA alignment tool. The same reference, GRCh37(used by 1000 genome) and human reference mtDNA sequence rCRS (revised Cambridge reference sequence), were used for ATAC-seq data processing. Samtools was used for manipulating sequence reads and calculating sequence depth. For all the data sets, the aligned reads were further filtered with mapping quality (Q > 30) and PCR redundancy was removed. The percentage of reads from mitochondrial genome compared to that of the nuclear genome were calculated after all the clean-up steps. The mitochondrial genome coverage was calculated using bases with sufficient sequence quality score (q > 30). A strong depletion region around 3107 due to the sequencing error(3170N) in the reference genome was excluded in the coverage plot.

23Most of the ATAC-seq pipelines remove mtDNA during their process. To rescue the genetic information from mtDNA, we modified our ATAC-seq pipeline and added SNP calling steps, which focuses on the mitochondrial genome. Briefly, adaptor sequences were trimmed from FASTQs using custom Python scripts. Paired-end reads were aligned to the reference genome using BWA. To improve the accuracy of heteroplasmic mutation calling, we followed the somatic mutation calling guidelines from GATK, with additional clean-up steps before variant calling. Reads mapped to mtDNA were extracted using Samtools from the final bam files and variants were called using VarScan2 with ‘--min-var-freq 0.001’ (Figure 1—figure supplement 1A). The heteroplasmic variants were further filtered through the following steps to exclude potential sequencing or mapping errors:

241. Thirteen frequent false-positive variants by misalignment due to extensive level of homopolymers in rCRS and due to sequencing error in the reference genome (reported in the previous study), were also observed and removed in this study. The following sites were explicitly removed:

25Misalignment due to ACCCCCCCTCCCCC (rCRS 302–315)

26A302C, C309T, C311T, C312T, C313T, G316C

27Misalignment due to GCACACACACACC (rCRS 513–525)

28C514A, A515G, A523C, C524G

29Misalignment due to 3107N in rCRS (ACNTT, rCRS 3105–3109)

30C3106A, T3109C, C3110A

312. Strand imbalance is a potential feature of sequencing error with various causes. To remove the potential sequence error from Illumina NextSeq (with a known high error rate at A bases) and sequence error from DNA damage(G->T, C->A), we required >2 reads detected from both the forward and reverse orientation, and strand is balanced (30% < forward/(forward +reverse)<70%).

323. Variant sites with VAF >0.9, but less than 1, were counted as homoplasmic variants. Although the germline polymorphic can be a backward heteroplasmic mutations, the observation of these events is higher than expected, which implies the false positive calling due to mapping bias for non-reference allele and sequencing errors.

334. For bulk ATAC-seq data from AML patients, heteroplasmic mutations with variant allele frequency >1% were reported.

34For all the AML cases (n = 15) from, we selected the cases (n = 12) with at least one confident heteroplasmic mtDNA mutation detected in any cell type for lineage relationship comparison. We found that in one patient (SU209), the number of heteroplasmic mutations (37) and their VAF are significantly higher than other patients. Most of these heteroplasmic mutations also overlapped with common variants present in the general human population (http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chrMT.phase3_callmom-v0_4.20130502.genotypes.vcf.gz), which indicates potential sample contamination. Therefore, this case was excluded from lineage relationship comparison and 11 AML cases were finally shown in Figure 2—figure supplement 1.

35To better evaluate the detection rate in single cell ATAC-seq data, we re-sequenced the previous libraries (LSCs and AML blasts from SU070 and SU373) from. The re-sequenced data were uploaded to GEO and accession number is GSE122576.

36Human AML samples were obtained from patients at the Stanford Medical Center with informed consent, according to institutional review board (IRB)-approved protocols (Stanford IRB, 18329 and 6453). Mononuclear cells from each sample were isolated by Ficoll separation, resuspended in 90% FBS +10% DMSO, and cryopreserved in liquid nitrogen. All analyses conducted here on AML cells used freshly thawed cells.

37Cell samples were first thawed and incubated at 37°C with 200 U/mL DNase in IMDM +10% FBS. To enrich for CD34 +cells, magnetic bead separation was performed using MACS beads (Miltenyi Biotech) according to the manufacturer’s protocol.

38For cell staining and sorting, the following antibody cocktail was used with the schema shown in Figure 2—figure supplement 3:

39CD34-APC, clone 581, Biolegend, at 1:50 dilution.

40CD38-PE-Cy7, clone HB7, Biolegend, 1:25 dilution.

41CD19-PE-Cy5, clone H1B9, BD Biosciences, 1:50 dilution

42CD20-PE-Cy5, clone 2H7, BD Biosciences, 1:50 dilution

43CD3-APC-Cy7, clone SK7, BD Biosciences, 1:25 dilution

44CD99-FITC, clone TU12, BD Biosciences, 1:20 dilution

45TIM3-PE, clone 344823, R and D Systems, 1:20 dilution

46CD45-KromeOrange, clone J.33, Beckman Coulter at 1:25 dilution

47Samples were sorted using a Becton Dickinson FACS Aria II. pHSCs were re-suspended and kept in cold FACS buffer containing 1 ug/mL propidium iodide prior to and after sorting. Cells were then immediately prepared for single cell ATAC-seq.

48Cells were washed two times in C1 DNA Seq Cell Wash Buffer (Fluidigm).~10 K cells were then re-suspended in 6 mL of C1 DNA Seq Cell Wash Buffer, and were combined with 4 mL of C1 Cell Suspension Reagent, 7 mL of this cell mix was loaded onto the Fluidigm IFC. Cells at a concentration of 260–380 cells/µL were then assayed using scATAC-seq as previously described. Briefly, single cells were captured using the C1 Single-Cell Auto Prep IFC microfluidic chips. Cells were permeabilized and accessible fragments were captured using 20 µL of Tn5 transposition mix (1.5x TD buffer, 1.5 µL transposease (Nextera DNA Sample Prep Kit, Illumina), 1x C1 Loading Reagent with low salt (Fluidigm), and 0.15% NP40) at 30 min at 37°C. In a 96-well plate, 7 µL of harvested libraries were amplified in 50 µL PCR for an additional 17 cycles (1.25 µM custom Nextera dual-index PCR primers in 1x NEBnext High-Fidelity PCR Master Mix using the following PCR conditions: 72°C for 5 min; 98°C for 30 s;) using the following PCR conditions: 72°C for 5 min; 98°C for 30 s; and thermocycling at 98°C for 10 s, 72°C for 30 s, and 72°C for 1 min. The PCR products were pooled creating a final volume of ~4.8 mL. The pooled library was purified on a single MinElute PCR purification column (Qiagen). Libraries were quantified using qPCR prior to sequencing. The scATAC-seq libraries were sequenced by Illumina MiSeq. The sequence data were uploaded to GEO under the accession number GSE122577.

49Single cell ATAC-seq were processed similarly to the bulk ATAC-seq, taking each individual cell as one sample. Recalibration steps were not applied for single cell data, as the sequence depth is not sufficient to empirically adjust the quality scores. After cleaning the alignment, files from every single cell were merged and heteroplasmic variants were first called with the merged bam and filtered using the same criteria as bulk data. Heteroplasmic variants called from merged data or from bulk data were re-counted in each individual cell using Samtools with ‘-q 20 -Q 20’. And the non-reference allele had to match the variants detected in merged or bulk data.

50In every single cell, if the variant allele detected in merged or bulk data were supported by any reads, it was considered positive; otherwise, it was counted as zero. A binary matrix was used to present the lineage relationship among single cells and plotted as a heat map. The intersections of the variants were quantified by the Upset R package. The number of detected variants showed a correlation with sequencing depth and the number of cells with all variants (Figure 1—figure supplements 3 and 4) confirmed the variants already co-existed in the ancestral cell. Following this assumption, the detection rate can be measured as the proportion of cells with variants in the total number of cells. For each variant, cells were separated into different bins, increased by 10, according to the total sequencing depth at each variant. The detection rate for each variant site was then calculated in each bin. The combined detection rate was estimated by 1-(1 R1)*(1 R2)*(1 R3)*(1 R4), where Rn is the detection rate for each variant.

51The probability of observing a mutation at a given site is Pn = n*r, where r is the average mutation rate in the mitochondrial genome and n is the copy of mtDNAs in a single cell. r is estimated to be ~10^-7 per base, n is around 100 ~ 10000 per cell, so Pn will be 10^-5 ~ 10^-3. The probability of N cells sharing the same mtDNA mutations, but arising independently, will be (Pn)^N. Thus, when there are more than 3 cells in the population sharing a common mtDNA mutation, the probability of these independently occurring will be close to 0. Cells with common mtDNA mutations having inherited the mutations from the same ancestral cell is more likely to explain the observed result. Furthermore, when a set of mutations (more than 1) is detected in more than 1 cells, the null hypothesis (independent occurrence) is rejected more confidently. The mutations within the ancestral cells can be inferred from the intersection of mutations. If a set of mutations co-existed in the ancestral cell, the absence of one of several linked mutations in the daughter cells is more likely caused by false negative detection in single-cell libraries or genetic draft during cell replication. Then the observed cells with different intersections (e.g V1 +V2) are expected to be given by Pv1*Pv2*N, after normalization by sequencing depth. The divergence of intersections from high-frequency mutations indicates the separation of mtDNA mutations and multiple cell lineage. The intersections of the variants were quantified by the Upset R package. In the scATAC-seq from pHSCs from SU353, the intersection of variants showed most of the cells were separated by two sets of different variants (Figure 2D). But there are a few cells displaying a mixture of variants from the two sets. We suspected these may cause by the doublet of cells in the same well during single cell separation on C1 chip. We further separated the intersection map by the chip and observed the number of cells with mixture variants correlated to the concentration of cells loaded to C1 Chip. These cells were removed during subsequent analysis. Single cells with any variants in the two sets were kept and cells with more than 40X coverage of mtDNA, but with no variants in the two sets were considered as HSCs with WT mtDNA. After all the filter steps, 153 pHSC cells had lineage information and were separated into three subgroups.

52ATAC sequences mapped to the nuclear genome were used for chromatin accessibility profiling. Bam files were merged for the same cell types and used as input files for chromVAR. Peak files from were used as open background regions to quantify the accessibility signal from every single cell. Cells with fewer than 200 unique reads or less than 25% of reads in peak regions were removed for chromatin analysis. chromVAR was applied to calculate TF motif-associated chromatin accessibility landscape changes and identify potential regulators of epigenomic variability. This approach quantifies accessibility variation across single-cells by aggregating accessible regions containing a specific TF motif, then compares the observed accessibility of all peaks containing a TF motif to a background set of peaks normalizing for known technical confounders. For determining differentially accessible motifs between different subpopulations, a Wilcoxon test was used to calculate the p values of the difference between the two groups.

53Custom analysis code can be downloaded from GitHub (https://github.com/ChangLab/ATAC_mito_sc;; copy archived at https://github.com/elifesciences-publications/ATAC_mito_sc).