Samples, mtDNA Enrichment, and Sequencing.

We studied the prevalence of mtDNA heteroplasmy in blood and buccal cells from 39 mother–child pairs residing in central Pennsylvania, analyzing 156 samples (39 mothers × 2 tissues + 39 children × 2 tissues) grouped in sets of four (two tissues from a mother and two tissues from her child). Total genomic DNA was isolated from each sample. Haplogroup analysis conducted via Sanger sequencing of the D-loop indicated European ancestry for all families (Dataset S1, Table S1). For each sample, we amplified mtDNA from total DNA in two overlapping 9-kb fragments and sequenced them with paired-end 250-bp reads on a MiSeq instrument (Materials and Methods). This enriched for mtDNA and minimized the presence of numts (Fig. S1 and SI Materials and Methods), the majority of which are short (15). Multiplexing 12 samples per run resulted in ∼10⁶ read pairs per sample. We confirmed the efficacy of our approach for mtDNA enrichment by applying it to Rho0 cells not harboring mitochondria (Fig. S2). To minimize potential contamination among samples, we followed previously devised guidelines (16) and used pUC18 and PhiX174 as spike-ins (Materials and Methods).

Heteroplasmy Discovery.

The sequencing read pairs were mapped to human mtDNA and nuclear genomes. On average, 85% of the reads per sample mapped to mtDNA (97% for samples without a spike-in; Fig. S3). We then applied our heteroplasmy discovery pipeline (Materials and Methods and Fig. S4). We required both reads of the pair to map uniquely, and in a proper orientation, to the reference mtDNA (Fig. S5). To compute minor allele frequency (MAF) at each site for each sample, we used bases with sequencing quality ≥30 from reads with mapping quality ≥20 (other thresholds led to almost identical results; Fig. S6). The mean sequencing depth per sample (averaged across sites) was 19,789× ± 770× (mean ± SEM). Tabulating depth on a per-site basis, 90% of sites in the mitochondrial genome were sequenced at ≥7,858× per sample (Fig. S7). The proportion of spike-in reads aligning to their respective references was as anticipated, suggesting absence of contamination among adjacent samples (Fig. S3).

In the search for heteroplasmies, we first identified sites with MAF ≥1% in individual samples. The sequencing depth per site required to detect true heteroplasmies with MAF ≥1% over the base quality error (0.1% for Phred score 30) with 99% power is 839× per site (one-sided power calculation for one-sample proportion test). Conservatively, we rounded up the depth requirement to 1,000×. A detection limit of MAF ≥1% allows detection of inherited and de novo variants that pass through the bottleneck if its size is <100. Mutations with lower frequency are accounted for with population genetics modeling (discussed below). After filtering for potential sequencing artifacts (Dataset S1, Table S2 and Materials and Methods), we identified 174 point heteroplasmies distributed among 100 quartets—groups of site-specific heteroplasmy frequencies from two tissues of a mother and two tissues from her child (Dataset S1, Table S3). These heteroplasmies were found in 31 families (eight families had no heteroplasmies).

Statistical Validation of Point Heteroplasmies.

To validate heteroplasmies, we used a novel statistical method that identifies heteroplasmic sites via a likelihood function accounting for instrument sequencing and mapping errors (SI Materials and Methods). With this method, all 174 point heteroplasmies were significant (P < 0.0003 for each site; Dataset S1, Table S4). Additionally, the allele counts for all 174 heteroplasmies tested were significant (P < 0.0003 for 172 sites, and P < 0.03 for the remaining two sites; Table S4) based on the variability observed for the same position among all samples (17).

Experimental Validation of Point Heteroplasmies.

We used Sanger sequencing to test all point heteroplasmies with MiSeq MAF ≥10% (Sanger method detection limit, Fig. S8A and Dataset S1, Table S5) in at least one sample per family and the corresponding sites from the other samples from the same family (we always sequenced newly amplified fragments). In total, we examined 21 sites × 4 samples = 84 sites, 44 of which had MiSeq MAF ≥10% (Dataset S1, Table S6). The presence of heteroplasmy was successfully validated in all these 44 cases. Thus, our false-positive rate for detecting heteroplasmies with MAF ≥10% is below 0.023 (1/44). The MAFs from the MiSeq and Sanger methods were well correlated (R² = 75%; Fig. S9A).

A set of point heteroplasmies with MiSeq MAF <10% was analyzed with droplet digital PCR (ddPCR) (18), which can detect heteroplasmies with MAF >0.2% (Fig. S8 B and C and Dataset S1, Table S7). Here we analyzed point heteroplasmies with MiSeq MAF between 1% and 10% in at least one sample per family and the corresponding sites from the other samples of the same family, a total of 10 sites × 4 samples = 40 sites, 18 of which had MiSeq MAF ≥1% (Fig. S9B and Dataset S1, Table S8). When we assayed the original amplicons used for MiSeq sequencing, the presence of heteroplasmy was confirmed in all these 18 instances. However, when we reamplified mtDNA from these 18 samples, in two instances (site 11,616 in M203C5-ch and site 11,825 in M210-bl) ddPCR did not confirm the presence of heteroplasmy. Repeating amplification and ddPCR for a third time again did not detect heteroplasmy (Dataset S1, Table S8), suggesting PCR errors in the amplicons sequenced with MiSeq. Thus, our false-positive rate for detecting heteroplasmies with MAF between 1% and 10% is 0.11 (2/18). Overall, the MAFs from the MiSeq and ddPCR methods were well correlated for the sequenced and newly amplified amplicons (R² = 95% and 79%, respectively; Fig. S9B).

Distribution of Point Heteroplasmies.

After removing two sites that failed to validate with ddPCR (discussed above), we retained 172 point heteroplasmies in 98 quartets (Dataset S1, Table S3). We assumed that these 98 point mutations arose independently in the families analyzed (or in their maternal ancestors). Point heteroplasmies were found at 87 unique mtDNA positions (Fig. S10). Six positions (185, 189, 214, 215, 16,093, and 16,183) were heteroplasmic in multiple families (four, three, three, two, three, and two families, respectively; Dataset S1, Table S3), likely owing to high mutation rate at the D-loop (13). Each mother on average carried 1.13 ± 0.04 heteroplasmies in her blood. This value was similar for maternal buccal tissue and for buccal and blood tissues of children (Fig. S11). Among the 98 point heteroplasmies 96 were transversions, resulting in a transition-to-transversion ratio of 48 (Fig. S10 and Dataset S1, Table S3).

There were significantly more and significantly fewer heteroplasmies in the D-loop and protein-coding regions, respectively, than expected based on their length and assuming equal propensity to harbor a heteroplasmy along mtDNA (Table 1). A high mutation rate for the D-loop had been documented previously (13). The nonsynonymous-to-synonymous rate ratio (d_N/d_S) at (concatenated) protein-coding genes was significantly lower than 1 (P = 5 × 10⁻³; Fisher’s exact test; Dataset S1, Table S9), suggesting purifying selection (19). Most nonsynonymous mutations were predicted to affect protein function (Dataset S1, Table S10).

Disease-Associated Mutations and Mutation Burden.

Eight families harbored eight point heteroplasmies (one per family) that can cause disease when present at high allele frequencies (Table 2). Among 39 mothers, 5 (or 1 in 8) were carriers of disease-associated mtDNA mutations in at least one of the two tissues analyzed. Mutations at four of the eight sites are associated with disease when homoplasmic for the mutant allele (20⇓⇓–23); however, in our data these were heteroplasmic (Table 2). For the other four of the eight sites above, disease can develop even when mutant alleles are heteroplasmic—with disease severity depending on the allele frequency. For A1555G, G13708A, and G3242A mutations, the allele frequencies were much lower than disease-associated frequencies (Table 2) (24⇓–26), suggesting lack of symptoms. Mutations at tRNA-Leu sites 3,242 and 3,243 contribute to several mitochondrial diseases (25, 27); notably, allele frequencies observed at site 3243 in the child of family M512 (Table 2) were comparable to those observed in mitochondrial encephalomyopathy, lactic acidosis, and stroke-like episodes patients (27).

Transmission of Heteroplasmies.

Considering heteroplasmy-containing quartets (Dataset S1, Table S3), we used the presence of heteroplasmy with MAF ≥1% in at least one sample from a family as a “prior” to support the existence of heteroplasmy at the same position for other samples of the same family if their MAF was ≥0.2% (greater than or equal to twice the value of the allowed sequencing quality error of 0.1%—Phred score 30). We classified 98 quartets into five categories (Table 3 and Dataset S1, Table S3) based on whether heteroplasmies were present in (i) both tissues of a mother and both tissues of her child (category “all,” which included dramatic shifts in allele frequency from mother to child, suggesting the germ-line bottleneck); (ii) both tissues of a mother, but absent from both tissues of her child (category “mother,” suggesting loss of a variant in the child owing to the germ-line bottleneck); (iii) both tissues of a child, but absent from both tissues of a mother (category “child” with candidate germ-line de novo mutations); (iv) both tissues of a mother and one tissue of a child, or in one tissue of a mother and both tissues of a child (category “somatic loss,” suggestive of a change in MAF in tissues owing to mitotic segregation) (12); and (v) one tissue of one individual of a family (category “somatic gain” with candidate somatic de novo mutations).

Site 4191 in family M500 seemed to harbor a de novo mutation in the child. ddPCR confirmed complete absence of the mutant allele in both maternal tissues but presence in both tissues of the child (with MAF of 4.7% and 5.6% in buccal and blood tissues, respectively; Table 3 and Dataset S1, Table S8). Examination of hair from the same individuals indicated homoplasmy in the mother and MAF of 1.2% in the child (Dataset S1, Table S11), confirming emergence of a novel allele.

The changes in allele frequencies between tissues of an individual, or between two generations, tabulated for our 98 quartets (Dataset S1, Table S3) followed an approximately normal distribution with mean zero (Fig. S12), corroborating the action of genetic drift as the major force affecting heteroplasmy allele frequencies (28). A decrease in allele frequency for a variant from mother to child will be indicative of purifying selection (29). When we plotted the relative change in allele frequency between mothers and children (Fig. S13), such a decrease was significant for nonsynonymous sites (P = 9.54 × 10⁻⁷, one-tailed nonparametric sign test), suggestive of purifying selection. Consistent with selection operating against transmission of nonsynonymous mutations, we observed a significantly lower proportion of these mutations among transmitted heteroplasmies (5 out of 43, or 12%, in “all” and “somatic loss” categories) compared with untransmitted heteroplasmies (8 out of 22, or 36%, in “mother” category; P = 0.025, Fisher’s exact test).

Comparing MAFs Between Tissues and Generations.

Compared with mtDNA in maternal tissues, mtDNA in child tissues underwent fewer mitotic segregations and replications and was exposed to mutagens for a shorter time. Therefore, we expect heteroplasmy allele frequency at a site to diverge less in the tissues of a child than in those of a mother. Indeed, the allele frequencies for the sites tabulated as quartets (Dataset S1, Table S3) were more strongly correlated between the two tissues for children (R² = 92%, Fig. 1A) than between the two tissues for mothers (R² = 49%, Fig. 1B). Stronger correlation for allele frequencies was observed between two tissues of a mother or of a child (discussed above) than between a mother and a child for the same tissue (R² = 13% for buccal, Fig. 1C; R² = 29% for blood, Fig. 1D), likely owing to the stronger action of the mtDNA germ-line bottleneck relative to mitotic segregation.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Correlation in heteroplasmy allele frequencies between (A) the two tissues of children, (B) the two maternal tissues, (C) buccal tissues of mothers and children, and (D) blood of mothers and children.

Maternal Age Effect.

We explored the relationship between age and the total number of point heteroplasmies for each individual. No association was found for children. For mothers we found a significant positive association (P = 0.039, 0.049, and 0.055 for combined, buccal, and blood heteroplasmies, respectively, Poisson regression; Fig. 2 and Fig. S14). Thus, older mothers accumulate more mutations in their somatic tissues, with the number of point heteroplasmies tripling over 30 y of life. Intriguingly, a positive association exists between the number of heteroplasmies in children and maternal age at fertilization (P = 0.010, 0.005, and 0.006, for combined, buccal, and blood heteroplasmies, respectively; Fig. 2 and Fig. S14). This suggests that older mothers accumulate more mutations in their germ-line tissues. In our dataset, there was a correlation between maternal age at fertilization and maternal age at sampling (R² = 43%, P = 0.002, linear regression), and we found that older mothers, who also had children later, likely transferred a larger number of accumulated mutations to their children (Fig. 2); whereas mothers who conceived under the age of 20 transmitted zero to one heteroplasmies, this number was two to three for mothers conceiving in their late 30s.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Maternal age effect. The dependence of the total number of heteroplasmies (blue) in mothers on their age at collection and (red) in children on maternal age at fertilization. Poisson generalized linear model fitted curves are indicated.

Estimating the Size of the Germ-Line mtDNA Bottleneck and Mutation Rate.

Because most heteroplasmy allele frequency changes between the two generations are consistent with genetic drift (Fig. S12), we can estimate the effective size of the germ-line bottleneck, that is, the size of the bottleneck in a traditional population model required to explain the observed amount of genetic drift [the actual number of mtDNA molecules passing through the bottleneck might be different, because they might segregate in units (8)]. Following the method developed by Millar, Hendy, and coworkers (30, 31), we assume that a child samples mutant mtDNA alleles at a given site from a binomial distribution with parameters p, the maternal MAF in the germ line (estimated here from somatic tissues), and N, the germ-line bottleneck size (SI Materials and Methods). Then the variance of the child’s heteroplasmy frequency at conception, or genetic variance, is σ²_gen = p(1 − p)/N. Solving for N, we obtain N = p(1 − p)/σ²_gen. We estimate the genetic variance as σ²_gen = σ²_raw − 4σ²_measure, where σ²_raw is the squared difference between the maternal and the child MAF at the site and σ²_measure is the uncertainty in measuring heteroplasmy frequency (includes sampling, PCR, and sequencing errors), which we estimated from sequencing amplified D-loop–containing clones (SI Materials and Methods; σ²_measure was multiplied by 4 because we are taking four measurements). This procedure produced an estimate of N in a quartet. To minimize false positives, we applied this approach to quartets where heteroplasmy was present in both maternal tissues (51 quartets in which at least one tissue in the mother had heteroplasmy with MAF ≥1% and the other tissue had MAF ≥0.2%; Dataset S1, Table S3), and thus likely was present in the maternal germ line. The median estimated N across these 51 quartets, when MAFs were averaged between the two maternal tissues and (separately) between the two child tissues, was 32.3 [interquartile range (IQR) 10.5–103.3; Fig. S15]. Similar estimates of bottleneck size were obtained when only blood (median N = 33.5, IQR 14.1–79.6) or only buccal tissues were used (median N = 29.8, IQR 9.5–68.1), and when quartets with nonsynonymous mutations were excluded (median N = 31.9, IQR 8.8–99.2). Accounting for the variance owing to mitotic segregation (SI Materials and Methods) led to median N = 35.0 (IQR 10.0–141.4; Fig. S15). Also, assuming that the single germ-line mutation we observe (at site 4191 in family M500, Table 3) originated in a single mtDNA segregating unit in the maternal germ line, and that its MAF in the child’s zygote was 3.8% (averaged across three tissues), we can estimate N as 1/0.038 = 26.3.

Next, we estimated the mtDNA germ-line mutation rate μ as in Millar, Hendy, and coworkers (30, 31). Assuming that new mutations enter the germ line at rate α, and that they are neutral and have equal probability to be transmitted to the next generation, only 1/N of them will go to fixation, leading to μ = α/N. A heteroplasmy can only be observed when its MAF is above a detection threshold θ. Analytically, it was shown (30, 31) that most heteroplasmies are lost without reaching θ, that most heteroplasmies reaching θ do not go to fixation, and that the rate of observed heteroplasmies can be approximated as μ₀ = 2αln(1/θ − 1). Solving for α, one obtains α = μ₀/(2 ln(1/θ − 1)). Thus, μ = μ₀/(2N ln(1/θ − 1)). Setting θ = 0.01 (our detection threshold) results in μ = 0.109 μ₀/N. Having observed 51 germ-line point heteroplasmies among 39 mothers, we estimated μ₀ as 51/(39 × 16,569 bp) = 7.9 × 10⁻⁵ heteroplasmies per transmission per site. Using N = 32.3, we thus estimated the mutation rate μ as 2.7 × 10⁻⁷ mutations per site per generation (IQR 8.3 × 10⁻⁸ to 8.2 × 10⁻⁷), or, assuming a generation time of 20 y, 1.3 × 10⁻⁸ mutations per site per year (IQR 4.2 × 10⁻⁹ to 4.1 × 10⁻⁸). The mutation rate estimate excluding nonsynonymous sites was 4.4 × 10⁻⁷ mutations per site per generation (IQR 1.4 × 10⁻⁷ to 1.6 × 10⁻⁶), or 2.2 × 10⁻⁸ mutations per site per year (IQR 7.0 × 10⁻⁹ to 7.9 × 10⁻⁸). That for the D-loop was 1.5 × 10⁻⁶ mutations per site per generation (IQR 4.8 × 10⁻⁷ to 4.7 × 10⁻⁶), or 7.7 × 10⁻⁸ mutations per site per year (IQR 2.4 × 10⁻⁸ to 2.4 × 10⁻⁷).

Indel Heteroplasmies.

Using the same thresholds (MAF ≥1% and depth >1,000×), we identified 120 instances of small indels affecting 10 unique mtDNA sites and 28 families (Dataset S1, Table S12 and SI Materials and Methods). All indels occurred in repeats—eight in homopolymer runs, one in a 9-bp tandem repeat, and one in a CA repeat. The latter indel was validated with ddPCR (Dataset S1, Table S13). The MAFs of indels in our samples were above microsatellite sequencing errors for our long-range PCR protocol (SI Materials and Methods). Further experimental validation will allow us to determine indel MAFs more accurately.