[1] 11w Detecting Heteroplasmy from High-Throughput Sequencing of Complete Human Mitochondrial DNA Genomes

[1] 251w Heteroplasmy, the existence of multiple mtDNA types within an individual, has been previously detected by using mostly indirect methods and focusing largely on just the hypervariable segments of the control region. Next-generation sequencing technologies should enable studies of heteroplasmy across the entire mtDNA genome at much higher resolution, because many independent reads are generated for each position. However, the higher error rate associated with these technologies must be taken into consideration to avoid false detection of heteroplasmy. We used simulations and phiX174 sequence data to design criteria for accurate detection of heteroplasmy with the Illumina Genome Analyzer platform, and we used artificial mixtures and replicate data to test and refine the criteria. We then applied these criteria to mtDNA sequence reads for 131 individuals from five Eurasian populations that had been generated via a parallel tagged approach. We identified 37 heteroplasmies at 10% frequency or higher at 34 sites in 32 individuals. The mutational spectrum does not differ between heteroplasmic mutations and polymorphisms in the same individuals, but the relative mutation rate at heteroplasmic mutations is significantly higher than that estimated for all mutable sites in the human mtDNA genome. Moreover, there is also a significant excess of nonsynonymous mutations observed among heteroplasmies, compared to polymorphism data from the same individuals. Both mutation-drift and negative selection influence the fate of heteroplasmies to determine the polymorphism spectrum in humans. With appropriate criteria for avoiding false positives due to sequencing errors, next-generation technologies can provide novel insights into genome-wide aspects of mtDNA heteroplasmy.

[1] 84w The mtDNA genome remains one of the most widely studied DNA segments in humans. It is particularly useful for studying population and evolutionary genetics because of its abundance in human cells, its uniparental, nonrecombining mode of inheritance, and its high mutation rate compared to that of the nuclear genome. 1 Although each individual is typically characterized by a single mtDNA type, in fact each individual is a population of mtDNA genomes, and the presence of multiple mtDNA types within an individual is termed heteroplasmy.

[2] 286w Although little noted at the time, the first report of heteroplasmy in humans was in 1983, involving a study of a noncoding region of human mtDNA from 11 human placentas. 2 Heteroplasmy has been investigated most often in correlation with mitochondrial disease, aging, and cancer. [3][4][5][6] To date, more than 400 mtDNA mutations have been associated with human disease, and most were observed in heteroplasmic states, with pathogenic mutations coexisting with normal mitochondrial genomes. 7 This suggests that the heteroplasmic level is of particular interest, as the disease phenotype becomes evident only when the percentage of mutant molecules exceeds a critical threshold value. Although this value differs for different mutations and in different tissues, it is usually in the range of 70%~90%. 8,9 Originally, heteroplasmy was believed to be quite rare in healthy individuals, 10,11 but subsequent studies found many non-disease-related heteroplasmies. [12][13][14][15] Moreover, heteroplasmy has started to play an important role in some forensic investigations. 16,17 Thus, heteroplasmy can also be a useful genetic marker. Regarding heteroplasmy as the intermediate stage between the generation of mutations and the fixation of mutations in the individual or cell, it represents polymorphisms within the populations of mitochondrial genomes in one cell or tissue. Thus, it can be a potential resource for studying the mutational pattern, possible role of natural selection, and even the existence of recombination in mtDNA. 18 For example, de novo mtDNA mutations in cancer tissues preferentially locate at the same positions as ancient variants in the human phylogeny, indicating similar selective constraints. 19 Understanding the basis, extent, and forces influencing the occurrence and subsequent fate of heteroplasmic mtDNA mutations is one of the principal challenges facing scientists and clinicians in the field of mitochondrial genetics.

[3] 202w A variety of techniques have been employed for heteroplasmy detection, including Sanger capillary sequencing, 13 high-performance liquid chromatography (HPLC), 20 pyrosequencing, 21,22 SnaPshot, 23 high-resolution melt (HRM) profiling, 24 a temporal temperature gradient gel electrophoresis (TTGE) strategy, 25 the Invader assay, 26 an amplification refractory mutation system, 27 and surveyor nuclease. 28 However, all of these methods have disadvantages, including the following: for some methods, the candidate heteroplasmic position needs to be defined first; the method may not allow determination of the actual heteroplasmic position; the level of heteroplasmy cannot be quantified accurately; and/or the method is too labor intensive to be applicable to large numbers of samples. Moreover, the efficiency of detection can vary substantially from laboratory to laboratory, even when the same method is applied, as a result of different instruments, chemistries, or standards for calling heteroplasmy. 29,30 In addition, most of the previous studies have restricted their examination to the control region or the hypervariable segments thereof; very few studies have analyzed heteroplasmies in the coding region, 31 even though these are likely to be of more importance for disease association. Therefore, more accurate and efficient methods are needed for the examination of heteroplasmy across the entire mtDNA genome.

[4] 102w Next-generation sequencing technologies can, in principle, provide the needed data, because millions of DNA reads can be produced in a single run at low cost. 32 This technology has been widely used in the detection of SNPs on a genome-wide scale, [33][34][35][36][37] and in one such study, 34 pooled DNA from 66 individuals was used to successfully detect SNPs. More interestingly, this study also found that the allele frequency inferred from the reads had a significant correlation with that inferred from genotyping. Such DNA pools resemble the mtDNA heteroplasmy scenario and therefore suggest that this approach should be useful for heteroplasmy detection.

[5] 139w However, the relatively high sequence error associated with next-generation sequencing technologies can produce false positives (i.e., sequencing errors may be falsely considered to be heteroplasmies). Moreover, problems could arise if PCR was involved during preparation of the sequencing library, because PCR may alter the proportion of different alleles and amplification error may also result in false positives. In this study, we used simulations, phiX174 control sequence data, artificial mixtures, and replicates to design and evaluate criteria for accurate detection of heteroplasmic positions with the Illumina GA platform. We then applied these criteria to a data set of complete mtDNA genome sequence reads for 131 individuals from five Eurasian populations that had been generated via a parallel tagged approach and sequenced with the Illumina GA platform. The results of this large-scale investigation provide insights into genome-wide patterns of mtDNA heteroplasmy.

[1] 90w Sequencing Overview A total of 45 million reads were generated for the 131 samples and four replicates, and 94% of the reads could be mapped to the rCRS. For the 36 bp reads, the average coverage was 67 5 7 (mean 5 standard deviation); meanwhile, 96.8% of the positions had coverage greater than 363 and at least ten reads for each strand. For the 76 bp reads, the average coverage was 211 5 21, and 98.6% positions had coverage greater than 363 and at least ten reads for each strand.

[2] 181w We used simulations to explore the effects of different sequencing error rates, coverage, and threshold of the minor allele frequency to define heteroplasmy, on the false-positive error (i.e., calling a sequencing error a heteroplasmy) and false-negative error (i.e., failing to call a true heteroplasmy) rates. Results are shown in Figure 1; 363 coverage is close to the lowest coverage that we have in our samples, whereas 763 coverage is close to the average coverage in our samples. The results illustrate the tradeoff that occurs between the false-positive and false-negative rates. When the heteroplasmy threshold is set relatively high (more than 10%), there are no false positives, even with a high sequencing error rate, but the false-negative rate is quite high, indicating that many true heteroplasmies will be missed. When the threshold is set relatively low (10% or less), the false-negative rate becomes very low, but now the false-positive rate becomes unacceptably high. Acceptable false-positive and false-negative rates can apparently be achieved only if the sequencing error rate can be reduced to 0.1% (Figure 1) or if the average coverage is increased.

[3] 314w The phiX174 genome is routinely sequenced to very high coverage in a separate lane in each Illumina GAII run, in order to provide baseline error rates for the base-calling software. Because in theory the phiX174 data should come from a single pure strain, there should not be any polymorphism, although five SNPs have been reported (Illumina technical support). Thus, any ''heteroplasmic'' position after assembly is presumably caused by sequencing error. We therefore investigated other ways to reduce false positives by applying various criteria to the phiX174 reads. Reads were randomly retrieved from the control lane to simulate a sequencing depth of 64-fold, and this procedure was repeated 100 times for each analysis. The number of heterozygous positions (defined as positions with at least one read with an alternative base) after assembly is shown in Table 1. Sequencing error is indeed an important issue, because in the absence of any QS filtering, 11% (583 out of 5386) of the positions were affected, and some of them can even reach a minor allele frequency between 20% ~30%. QS filtering can improve the situation, but even a fairly stringent QS (QS R 23) leaves many heterozygous positions (Table 1). However, double-strand validation (i.e., requiring at least one read from each strand) considerably reduces the impact of sequencing error; with no quality filter, double-strand validation reduces the average number of heterozygous positions from 583 to 18 (Table 2). And with a combination of a stringent quality filter and double-strand validation, only two heterozygous positions are detected on average. Note that this better performance of double-strand validation was not caused by any bias in the number of reads in one direction, because all positions have a coverage greater than 40 and at least ten reads in each direction. Instead, sequence errors tend to be context specific, and hence the majority of sequence errors are found on only one strand.

[4] 143w To investigate the performance of double-strand validation with various frequency thresholds for calling heterozygous positions, we carried out additional simulations using an average coverage of 66, to mimic our real data. Table 2 compares the results of single-versus double-strand validation for various thresholds (for other sequencing depths, see Table S2), requiring QS R 20 at the position in question (and QS R 15 for the 5 bp flanking sequence on each side). As expected, a lower frequency threshold resulted in a higher false-positive rate, but only one false positive was consistently detected with double-strand validation and a frequency threshold of 10% or more. This position (1301) is not included in the five SNPs previously reported, but it has been found to be polymorphic by another laboratory (T. Skelly, personal communication), so it is likely to be a novel SNP in the phiX174 genome.

[5] 97w On the basis of both the simulations and the phiX174 analyses, we decided on the following criteria for calling a heteroplasmic position: QS R 20 at the position in question, QS R 15 at the 5 bp flanking either side, a minor allele frequency of at least 10%, and the minor allele observed in at least one read in each direction. These criteria should result in no false positives, and we should have 100% power to detect heteroplasmies present at a frequency of 20% or more and 50% power to detect heteroplasmies at a frequency of 10%.

[6] 165w To test the sensitivity and specificity of these criteria, we applied them to three artificially mixed samples consisting of DNA derived from two individuals in different ratios (1:9, 1:3, 1:1) whose mtDNA genome sequences differ at 25 positions. All 25 positions were observed to be heteroplasmic, with inferred frequencies closely corresponding to those expected (Figure 2), except for six positions for the 1:9 mixture, which fell close to, but below, the 10% threshold. The Pearson correlation coefficient between the heteroplasmic level estimated from sequencing and the actual proportion is 0.944 (p < 0.001), in good agreement with a previous study of artificially mixed samples analyzed by high-throughput sequencing. 47 Four potential false positives were observed: position 7030AC (the first letter denotes the majority nucleotide and the second letter the minority nucleotide), with a minor allele frequency of 12% in the sample with a 1:1 ratio; position 13604GA, with minor allele frequencies of 10% and 12% in the 1:3 and 1:9 mixtures, respectively; and the A

[7] 174w homopolymer region beginning at 12418, where 11% of the reads have 7A and 89% have 8A in the 1:9 mixture. However, as shown below, position 13604GA is actually inferred to be heteroplasmic, with a minor allele frequency of 11%, in the sample (Ir43) that was used as the major component in the mixed samples. Moreover, the A homopolymer region beginning at 12418 was heteroplasmic in both of the samples, with a minor allele (7A) frequency of 12% and 10%, respectively. Therefore, only position 7030 qualifies as a true false positive. A more stringent requirement of QS R 23 at the heteroplasmic position (rather than QS R 20) would eliminate this false positive, but it would also remove an unacceptably high number of other reads for the true heteroplasmies. However, requiring two distinct reads from each strand would also eliminate this false positive while retaining all of the other inferred heteroplasmies. Therefore, we adjusted the criteria for heteroplasmy to require that the minor allele was observed in at least two independent reads from each strand.

[8] 119w Frequency Variation between Replicates Four samples were sequenced twice independently, allowing us to assess how inferred minor allele frequencies varied between replicates. The average frequency variation between replicate pairs was 0.0057 5 0.011, and there were 42 positions having a minor allele frequency difference greater than 0.1 (Figure 3). More stringent QS filtering (QS R 23 for the position in question and QS R 20 for the 5 bp flanking sequence on each side) reduced but did not eliminate the discrepancies; 19 positions still had a minor allele frequency difference greater than 0.1. Moreover, low coverage does not appear to play a significant role, because only five positions have a coverage lower than 10 in one of the replicates.

[9] 177w However, for 21 of the 42 positions with a minor allele frequency difference bigger than 0.1, the nucleotides involved were A and C (Table S3). Moreover, for many of these positions, the minor allele is observed on only one strand, indicating that the discrepancy between replicates is due to an increase in the minor allele on one strand only, consistent with sequencing errors. This same pattern was observed for another four positions in which A and T were involved (Table S3). The one exception to this pattern is position 385 for Az41, for which A and G are observed on both strands in both replicates. The frequency of the minor (G) allele was 20% versus 8.7% in the replicates; however, the difference was not significant (p ¼ 0.059, Fisher's exact test). It is likely that this is a real heteroplasmic position and that the frequency difference between replicates is due to stochastic effects. All of the remaining 16 positions were located in the homopolymer region, for which inaccurate alignment can explain the large minor allele frequency difference.

[10] 163w Application of the criteria developed above to the full data set of 131 individuals resulted in the detection of 37 point heteroplasmies among 32 individuals (Table 3), or 24% of the individuals studied. The average coverage for the heteroplasmic positions is 79 (range: 28-170), and for only four positions (3014GT, 3492AC [twice], 10208CT) was the minor allele observed in fewer than six reads. Five individuals possessed two heteroplasmies each; overall, the number of heteroplasmies per individual did not differ significantly from the expectation when 37 heteroplasmies are randomly assigned to 131 individuals (p > 0.05). There were also no significant differences in how heteroplasmies were distributed among populations or among haplogroups. However, three heteroplasmies (146CT, 3492AC, 16223TC) occurred at the same position in two different individuals who belonged to different populations and had different mtDNA haplogroups; thus, these are likely to represent independent heteroplasmies. This is significantly more than expected by chance if heteroplasmies occur randomly across the mtDNA genome (p < 0.001).

[11] 151w We also detected three indel heteroplasmies (outside homopolymer regions) among three individuals (Table 4). One of these was located in a tRNA gene, whereas the other two were in the control region. In addition, we observed another 112 indels in the following homopolymer and STR regions (with the variation observed indicated in parentheses): 66~71 (5C-6C), 303~309 (6C-9C), 514~523 (5CA-6CA), 12418~12425 (7A-8A), and 16184~16193 (9C-12C). Even if we only use the reads that span these homopolymer and STR regions, we still observed these length variations in at least two distinct reads per strand. However, slippage of the DNA polymerase during replication can also result in length variation, which is frequently observed in the homopolymer and STR region. 48,49 Therefore, the indel heteroplasmic regions would need additional validation, ideally by a method that does not require PCR amplification, and for this reason we do not include these regions in our analysis of heteroplasmy.

[12] 67w We selected nine heteroplasmic sites for independent validation via single-base extension assays. For each site, the predicted heteroplasmy state was observed in the corresponding individual (Figure S1), whereas control individuals lacking the heteroplasmy displayed only single peaks in the assays (data not shown). These results indicate that the heteroplasmies inferred by analysis of highthroughput sequencing reads are indeed true heteroplasmies, because they are validated by another method.

[13] 98w Thirteen of the 34 heteroplasmic positions are located in the control region, which is more than expected by chance if heteroplasmic positions are distributed randomly across the mtDNA genome (p < 0.001, chi-square test). The ratio of transitions to transversions at heteroplasmies is 8.25, which is not significantly different (p ¼ 0.616, Fisher's exact test) from the ratio of 8.09 for polymorphic positions in the same individuals. Moreover, mutational spectrum for heteroplasmies does not differ significantly from that for polymorphisms in the same populations, or from the mutational spectrum for polymorphisms reported in the Mitomap database (Figure 4).

[14] 130w For investigation of the correlation between the incidence of heteroplasmy and mutation rates, the RMR for each heteroplasmic position was obtained from a previous study. 45 We classified the heteroplasmic positions into three categories: entire genome, control region, and coding region; the comparisons of RMR between heteroplasmic positions with all positions and all mutable positions in each category are shown in Table 5. Heteroplasmic posi- tions have higher mutation rates in all categories, significantly so except when compared to all mutable positions in the coding region. Overall, heteroplasmies are occurring preferentially at positions with high mutation rates. This is particularly evident in the control region, in which six of the ten positions with the highest mutation rates (Soares et al. 45 ) were found to exhibit heteroplasmy in our study.

[15] 94w We detected 21 heteroplasmies at 20 positions in the protein-coding genes (3492AC is shared by two individuals). Ten of these are nonsynonymous mutations, ten are synonymous mutations, and one is located in the overlap region involving the ATP6 (MIM 516060) and ATP8 (MIM 516070). There is an excess of nonsynonymous mutations at the heteroplasmic positions when compared to polymorphism data from the same individuals (164 nonsynonymous and 385 synonymous mutations; p ¼ 0.051, Fisher's exact test) or to a previous study 50 (413 nonsynonymous and 1037 synonymous mutations; p ¼ 0.045, Fisher's exact test).

[16] 271w Five of the ten nonsynonymous heteroplasmies are located at positions with an assigned RMR, 45 and one of them is located at the position with the highest RMR among 824 nonsynonymous mutations (5460AG, RMR ¼ 36). Two other nonsynonymous heteroplasmies also have high RMR (15314GA, RMR ¼ 12; 11253CT, RMR ¼ 5), whereas the other two (7754AG, 8743GA) both have RMR ¼ 1. The heteroplasmic nonsynonymous sites tend to have higher mutation rates than previously reported nonsynonymous mutations in polymorphism data 45 (p ¼ 0.068, Mann-Whitney test). The other nonsynonymous heteroplasmies-3492AC (which occurred twice), 3532AG, 14561AG, and 13604AG-have never been observed as polymorphisms (based on sequences retrieved from Mitomap, 7 mtDB, 43 and GenBank, as well as 700 unpublished mtDNA genome sequences from our laboratory). The overall ratio of nonsynonymous heteroplasmies at sites without previously reported mutations, 5/10, is significantly higher than that observed in polymorphism data from the same individuals, 23/164 (p ¼ 0.01, Fisher's exact test). There is thus a dichotomous tendency for nonsynonymous heteroplasmies to occur either at rapidly evolving sites or at sites that have not been observed as polymorphisms. Eight of the ten synonymous mutations have an assigned RMR; 45 the average RMR for these eight positions is 3.0 5 2.6, which is higher than the overall Finally, two heteroplasmic positions (3014GT, 1552AG) are located in the stems of rRNA genes, 51,52 with the major alleles identical to the rCRS, and were not observed to be polymorphic in the database. One heteroplasmic position (10427AG) is located in the connection of two stems of the tRNA-ARG (MIM 590005), 53 and this position has an RMR of 2.0.

[17] 80w Four of the 34 (11.8%) heteroplasmic positions are reported as disease-associated in Mitomap 7 (Table 3), which is more than expected by chance (2.5% of all mtDNA positions are reported to be disease-associated in Mitomap; 7 p ¼ 0.010, permutation test). However, all of these positions have also been reported previously as polymorphisms; moreover, there is also an excess number of polymorphisms at disease-associated sites in the mtDNA genome sequences from the same individuals (6.1%; p < 0.001, permutation test).

[18] 128w It seems likely that disease-associated mutations with mild effects will be observed as polymorphisms in ''normal'' individuals more frequently than disease-associated mutations with strong effects. Overall, 22.4% of the disease-associated mutations in Mitomap are also present in normal populations and 77.4% are absent in normal populations. All of the four disease-associated heteroplasmic positions are observed as polymorphisms in normal populations, which is significantly more than expected (p ¼ 0.003, permutation test). Moreover, the average RMR for these four disease-associated heteroplasmic positions is 51, which is significantly higher than the average RMR of 2.72 for all mutable positions (p < 0.001, Mann-Whitney test) and significantly higher than the average RMR of 8.8 for all disease-associated positions that also present as polymorphisms in normal populations (p ¼ 0.007, Mann-Whitney test).

[19] 100w Although the high-throughput sequencing technologies have a higher per-base sequencing error rate compared with the traditional capillary sequencers, the new technologies still give highly accurate sequences because of the higher sequencing depth. 39,42,54 However, this holds only for assembling a consensus sequence or calling SNPs, in which using a ''majority rule'' eliminates the impact of sequencing error and produces the desired result. Trying to detect heteroplasmy presents different issues, because the minor allele at the heteroplasmic position may be indistinguishable from sequencing error. QS filtering can help eliminate sequencing errors but is insufficient to completely solve the problem (Table 1).

[20] 399w The most obvious way to distinguish between sequencing errors and heteroplasmy is to invoke a threshold. Sequencing errors happen at a specific rate (10 À2 ~10 À3 , based on our unpublished data and other studies 39 ), so frequency of an error observed at low coverage should decrease to this level with higher coverage, whereas the level of heteroplasmy should not change with the coverage. However, sequencing errors are reported to be position dependent and context dependent. 36,39,42 In fact, we have observed several sequencing error ''hot spots'' (including 257AC, 3492AC, 3511AC, 4774TA, 5290AT, 9801GT, 10306AC, 10792AC, 11090AC), defined here as meeting the following criteria: (1) R 10% of the individuals have the minor allele with frequency R 10% at this position, but (2) the minor allele for R 90% of the individuals cannot be validated by reads from both strands. As expected, most of the sequencing-error hotspots involve misidentification of A to C; thus, some errors can reach a high enough frequency to pass a simple frequencythreshold criterion and contribute to false positives. However, given that heteroplasmy is not strand dependent and that it is highly unlikely that an error-prone context would exist for the same position on both strands, a useful filter is double-strand validation: heteroplasmies are called only if both of alleles can be detected by reads from both strands. Double-strand validation alone significantly reduced the number of heterozygous positions detected in the phiX174 sequence data (Table 1). However, even double-strand validation is not sufficient to eliminate all sequencing errors: one position (7030AC) in the artificially mixed samples passed this filter but is a false positive. This can occur because sequencing error also occurs in the other direction by the usual background error rate (e.g., with 40 reads in the non-error-prone direction and an error rate of 0.3%, there is an 11% chance of observing at least one error). We therefore implemented the requirement of observing both alleles in at least two independent reads from each direction to validate the heteroplasmy, and the number of reads required on both strands should be increased with a higher sequencing depth. Although this requirement appears to be sufficient to eliminate virtually all false positives due to sequence errors, it does increase the probability of false negatives due to failure to observe a true minor allele in at least two reads from each direction, especially when coverage is low.

[21] 174w Another aspect of Illumina GA sequencing errors is that they tend to be position dependent: typically, the error rate at the end of the read is approximately 2-to 6-fold higher than at the beginning. One potential way to reduce the false-positive errors that are caused by position-dependent error might therefore be to require that heteroplasmies are validated by the beginning or middle of reads. Indeed, this requirement significantly reduced the number of heterozygous positions in the phiX174 reads (see details in Table S4). However, our data also suggest that the error at the ends of the reads is equivalent to that in the middle of the reads under the QS filtering that we applied (QS R 20 for the position in question and QS R 15 for the 5 bp flanking sequence on each side), which may reflect the fact that QS is correlated with the position in the read. We therefore did not include the requirement of confirmation by the beginning or middle part of the reads in order to call a heteroplasmy.

[22] 92w In addition to sequencing errors, inaccurate assembly of the reads can create artificial heteroplasmies, because the minor allele could come from reads originating from another part of the mtDNA genome. To evaluate whether heteroplasmy is caused by improper alignment, we removed all reads with multiple best hits on the final consensus sequence (i.e., mapping quality equals 0 when assembling the reads by using Burrows Wheeler alignment as implemented in the software BWA 55 ), and none of the inferred heteroplasmies disappeared. Therefore, incorrect assembly does not appear to be influencing our results.

[23] 159w Another question concerning the assembly is whether to use all reads or only unique reads (i.e., those with different start and/or endpoints) from high-throughput sequencing data. [56][57][58] The concern is that duplicate reads may represent copies from the same molecule rather than independent reads. We examined this question in the artificially mixed samples and found that using unique reads versus all of the reads does not influence the detection and level of heteroplasmy (results not shown). However, when we resampled the reads from the 1:3 mixture to investigate the influence of sequencing depth, using duplicate reads resulted in higher false-positive and false-negative error rates than using unique reads (Table S5). Apparently, with low coverage even a small number of duplicate reads can have a large impact on the inferred minor allele frequency. Because a large fraction (~27%) of the positions in our samples have a sequencing depth less than 1503, we used only the unique reads to infer heteroplasmies.

[1] 20w We first discuss issues related to the criteria for detecting heteroplasmy, and we then discuss the implications of our results.

[1] 179w The mtDNA sequence reads used here were generated in a previous study (unpublished data). In brief, the entire mtDNA genome was amplified from 147 individuals from five populations (Georgia, Armenia, Azerbaijan, Iran, and Turkey) in two overlapping products of about 9.7 and 7.3 Kb and sequenced on the Illumina GAII platform (GAII; San Diego, CA, USA) via a multiplex sequencing protocol for sequencing libraries; 38,39 details are described elsewhere (unpublished data). Out of 147 samples, the reads from 131 met the criteria for accurate detection of hetero-plasmy, developed below. Of these, 97 were sequenced once with single-end 36 bp reads to an average coverage of 653, 17 were sequenced twice with single-end 36 bp reads (because of insufficient coverage from the first lane) to an average of 783, and 17 were sequenced once with single-end 76 bp reads to an average coverage of 2113. In addition, four samples originally sequenced once with 36 bp reads were resequenced (with the use of new PCR products amplified from the original DNA templates) with single-end 76 bp reads for assessment of reproducibility.

[2] 155w The reads were assembled with the software MIA, 40 the revised Cambridge reference sequence (rCRS) used as the reference. 41 MIA is optimized for circular genome assembling and performs gapped assembling iteratively: after a consensus is called, the reads are realigned to the consensus and a new consensus is called; this process iterates until it converges on a single consensus sequence. For the 36 bp reads, maximally, three mismatches or two mismatches plus one gap were allowed to successfully map the read, whereas for the 76bp reads, five mismatches or four mismatches plus one gap were allowed. In addition, any read with more than two bases having a Phred-like quality score (QS) lower than 15 was removed (for 76 bp length reads, the threshold was five bases with QS < 10), and duplicate reads (reads mapping to the same position with same orientation on the reference) were removed, keeping the one with the highest QS.

[3] 86w The template DNAs from two individuals differing at 25 positions in their complete mtDNA genome sequences were mixed in different proportions (1:1, 1:3, and 1:9). DNA concentrations were measured by a NanoDrop ND-1000 spectrometer (NanoDrop Technologies, Wilmington, DE, USA), diluted and mixed in the desired proportions, and then used for long-rang PCR amplification and Illumina GAII sequencing with 76 bp reads. These artificially mixed samples were used for examination of the correlation between the heteroplasmy level (mixture proportion) and minor allele frequency estimated from the reads.

[4] 210w To evaluate the impact of sequencing error on the detection of heteroplasmy, we performed simulations for varying levels of sequencing error, heteroplasmy level, and sequencing depth. For each simulation, one biallelic heteroplasmic position was assigned to the genome randomly, with the heteroplasmy level (minor allele frequency) set to be 5%, 10%, 20%, 30%, or 40%. Reads with a specific length (36bp, 76bp) were generated randomly from the rCRS to reach an average coverage of 363 or 763. The sequencing error in the simulation was evenly distributed along the reads and along the mitochondrial genome, with the error rate set to be 1%, 0.5%, 0.3%, or 0.1%. Because Illumina GA uses the same laser to excite two pairs of nucleotides (A/C and G/T), the pairs produce similar emission spectra and are thus poorly distinguished by optical filters, which results in a higher proportion of errors involving A/C and G/T. 39,42 On the other hand, heteroplasmy is highly biased toward transitions (A/G, C/T), as with mtDNA substitutions. 13,14 We therefore specified that all sequencing errors would be A/C or G/T changes, whereas all heteroplasmy would involve A/G or C/T changes, so that we could readily distinguish sequencing errors from heteroplasmy in our simulations. For each setting, the simulation was repeated 100 times.

[5] 70w Polymorphic positions were retrieved from Mitomap, 7 mtDB, 43 4775 whole mitochondrial genomes from GenBank, 44 and 700 whole mitochondrial genomes from our laboratory (unpublished data). The relative mutation rate (RMR) for each polymorphic position was retrieved from Soares et al. 45 Disease-associated positions were retrieved also from Mitomap 7 (April 27, 2010 version); when defining a heteroplasmy as disease-associated, both the nucleotide position and the specific mutation were considered.

[6] 36w Haplogroup-defining positions were retrieved from Phylotree, 46 and an in-house script was used to collect all the haplogroupdefining positions and the corresponding mutation type, which were then used to identify heteroplasmies occurring at such haplogroup-defining positions.

[7] 156w To independently verify the heteroplasmic positions inferred from high-throughput sequencing, we performed single-base Figure 1. False-Negative Error Rate and False-Positive Error Rate in Detecting Heteroplasmy Inferred from Simulation False-negative error rate and false-positive error rate calculated under different error rates (1%, 0.5%, 0.3%, 0.1%), coverage (363, 763), heteroplasmy levels (minor allele frequencies of 10%, 20%, 30%, 40%), and frequency thresholds used to define heteroplasmy (40%, 30%, 20%, 10%). For each setting, the simulation was repeated 100 times. FP denotes the false-positive error rate, FN denotes the false-negative error rate. extension assays for nine heteroplasmic sites. Primers were designed to amplify shorter products around these sites for both control and heteroplasmic individuals with the use of original DNA as template (Table S1, available online). All single-base extension assays were carried out with the ABI Prism SnaPshot Kit (Applied Biosystems), with amplicons detected via capillary electrophoresis on an ABI Prism 3100 Genetic Analyzer according to the manufacturer's instructions.

[8] 134w PCR is another factor influencing the detection of heteroplasmy; biased amplification can alter the minor allele frequency, and sequence error introduced during PCR amplification may be falsely detected as heteroplasmies. 59,60 However, several observations indicated that PCR has had little or no effect on our analyses: nine inferred heteroplasmy were genotyped by single-base extension assays, and all nine were validated; the inferred heteroplasmy level in artificially mixed samples was close to the mixture ratio; and sequencing of four samples in duplicate showed no frequency difference greater than 10% that could be validated by the double-strand criterion. It should be noted that all PCR products used in the above studies were independently amplified from the original DNA template. Therefore, possible artifacts from the PCR process do not appear to be influencing our inferences concerning heteroplasmy.

[9] 78w Other Factors that Influence the Detection of Heteroplasmy Besides sequencing error, assembly error, and PCR artifacts, there are other processes that could potentially influence the identification of heteroplasmic positions by generating false positives. These include jumping PCR, 61 inadvertent sequencing of nuclear mitochondrial pseudogenes (numts), 48,62 contamination during the library preparation, and the ''stochastic effect'' produced by random sampling of sequencing reads. To what extent could these be responsible for the heteroplasmies that we detected in this study?

[10] 55w The first possibility is jumping PCR, which could result in the incorrect allocation of a read to an individual, because a parallel tagged approach was used in our study. 61 However, after the indexing PCR of seven cycles, no additional PCR was performed on the pooled samples, so jumping PCR should not be a problem.

[11] 79w Numts are another potential source of contamination during mtDNA PCR amplification, because there are 46 paralogous nuclear DNA fragments that represent the entire mitochondrial genome. 62 However, instead of shortrange PCR or capture-based methods, which may be more severely influenced by numts, we utilized long-range PCR. No secondary amplification products were observed, and moreover, 94% of the reads could be mapped to the mtDNA genome, indicating that the reads are indeed originating from authentic mtDNA and not from numts.

[12] 107w Inadvertent sample mixtures or contamination could also produce the appearance of heteroplasmy (indeed, this was the basis for the artificially mixed samples). To verify that such contamination does not explain inferred heteroplasmies, we checked whether the content of the heteroplasmies in one individual could be used to define two different haplogroups in Phylotree. 46 In fact, several individuals in the original study (unpublished data) did exhibit numerous putative heteroplasmic positions that could be explained by sample mixtures involving two different haplogroups; these individuals were excluded from further analysis. In the 131 individuals used in this study, sample mixtures involving any known haplogroups cannot explain the inferred heteroplasmies.

[13] 95w A final concern is the ''stochastic effect,'' i.e., how much variation can occur in sample replicates, and how this would influence the inference of heteroplasmy. In the four replicates, everything was repeated from the template DNA via the same methods; thus, any allele frequency variation for each position between these replicate pairs should be caused by the stochastic effect. Although the ''stochastic effect'' could result in frequency variation greater than 0.1, none could be validated by the doublestrand criteria, indicating that stochastic effects are unlikely to produce false positives under our criteria for detecting heteroplasmy.

[14] 17w Supplemental Data include five tables and one figure can be found with this article online at http://www.cell.com/AJHG/.

[1] 125w In each Illumina GA sequencing run, the phiX174 bacteriophage genome is sequenced on a separate lane to very high coverage (about 60,0003) and then used as a control for sequencing-error correction for that run. Theoretically, there should be no heteroplasmy in the phiX174 genome, given that it is from a single strain, although in fact Illumina technical support reports five SNPs in phiX174 genome sequences. Therefore, we applied our criteria for detecting heteroplasmy to the phiX174 data, to estimate the false-positive error rate associated with various criteria. To simulate different sequencing depths, we retrieved subsets of the reads from the phiX174 control, and for each depth (15-, 30-, 60-, 90-, 120-, 150-, 180-, 210-, 600-, and 1200-fold average coverage), we repeated the simulations 100 times.

[2] 104w As discussed above, the heteroplasmies detected by our criteria are unlikely to be caused by sequencing errors or other artifacts. We now discuss the insights provided by this large-scale examination of mtDNA genome-wide heteroplasmy. We identified 37 point heteroplasmies and three indel heteroplasmies among 33 of 131 individuals (25%), which is higher than the range of 3.8%-6% reported in previous studies that used other methods. 13,63 We attribute this higher rate to the increased sensitivity that high-throughput sequencing offers for detecting heteroplasmy and to the fact that we are analyzing heteroplasmy across the entire mtDNA genome, whereas previous studies focused on the control region.

[3] 97w The familiar transition bias in human mtDNA muations 64 was observed, in that 89% of the point heteroplasmies are transitions. Thirteen point heteroplasmies (35%) were found in the control region, which is more than expected, and all of them are located in hypermutable positions (RMR R 6). An association between heteroplasmy and hypermutable positions was found previously, 65 and overall, these results suggest that mutation is the major driving force behind heteroplasmy and that a mutation-drift process can explain how heteroplasmies arise, drift to high frequencies within an individual, and eventually become ''fixed'' as polymorphisms among individuals.

[4] 133w However, in the protein-coding region, mutations are located in both hypermutable positions and hypomutable positions. Altogether, eight heteroplasmies are located in positions that have never been reported to mutate, which is more than expected by chance (p ¼ 0.04). Moreover, there is a significant excess of heteroplasmies involving nonsynonymous changes in comparison to polymorphism data. It thus seems as if deleterious mutations arise as heteroplasmies and can reach appreciable frequencies (> 10%) within individuals but do not drift to fixation. Instead, purifying selection must be operating on some heteroplasmies to prevent their fixation within individuals. Although purifying selection on mtDNA has been inferred from other studies, [66][67][68] no previous evidence of such selection involving heteroplasmy has been found, 69 which may reflect the limited number of heteroplasmies previously studied in the coding region.

[5] 132w Significantly more heteroplasmies were detected in disease-associated positions than expected by chance, but all of them were also observed to be polymorphic in the same population or in our polymorphism data set and, moreover, are associated with high mutation rates. These are not characteristics expected to be associated with deleterious mutations and hence may reflect limitations on the accuracy of studies that attempt to elucidate disease-associated mtDNA mutations. Or, these may be associated with mild effects on the phenotype. Still, we found that even in comparison to the polymorphism profile, there is an excess number of heteroplasmies in disease-associated positions. It is thus possible that these disease-associated heteroplasmies may drift to high enough frequencies within an individual to result in the disease phenotype and/or that they will be removed by purifying selection.

[6] 286w We have developed a set of criteria to detect heteroplasmy in complete human mtDNA genomes from the reads generated by Illumina GAII technology. Recently, another study used the same sequencing technology to investigate mtDNA heteroplasmy. 14 This study sequenced only one individual per lane, to an average coverage of ~16,000, and hence used a correspondingly lower threshold for the minor allele frequency of 1.6% to call heteroplasmic positions. They detected 40 heteroplasmies above this threshold in ten healthy individuals; 14 of these, seven heteroplasmies in four individuals had a minor allele frequency greater than the 10% threshold that we used. This is significantly more than the 37 point heteroplasmies that we detected among 131 individuals (p ¼ 0.02, resampling test). Possible reasons for this difference include: a higher false-negative error rate in our study, which is due to lower sequencing depth; the different tissues used (blood and saliva in our study, colonic mucosae in their study 14 ); and/or the relatively old age of individuals in their study. 14 Although sequencing to a much higher depth obviously increases the power to detect heteroplasmies by lowering the minor allele frequency threshold, stringent quality control procedures are still necessary to determine a proper frequency threshold. In fact, we note that in the previous analysis of heteroplasmic variants among tissues from a single individual (Table 2 of He et al. 14 ), five of 14 heteroplasmies (36%) involve AC/GT changes, which is the most common sequencing error on the Illumina GA platform. This is significantly more than the 9% frequency of AC/GT changes among heteroplasmies in our study (p < 0.05), and it therefore suggests that some of these are likely to reflect false positives due to sequencing errors.

[7] 203w The strategy used to detect heteroplasmic variants can therefore vary, depending on the goal of the study (and available funding). Sequencing one sample per lane (or even one sample on multiple lanes) will lower the threshold for detecting heteroplasmic variants but will increase the cost and therefore limit the number of samples that can be studied. Parallel tagged approaches, which were used here, result in lower coverage and hence require a higher minor allele frequency to call a heteroplasmy, but they are a more cost-efficient approach for analyzing heteroplasmy in a large number of samples. Indeed, we predict that the next few years will see a huge increase in complete mtDNA genome sequences generated to an average coverage of ~503-1003 with next-generation sequencing platforms. These data will provide a rich resource for further investigation of heteroplasmy. Regardless of the strategy and average coverage obtained, our results indicate that accurate calling of heteroplasmic positions requires the analysis of control data and an appropriate statistical model to generate appropriate criteria. With such appropriate criteria for avoiding false positives due to sequencing errors, high-throughput sequencing platforms can provide a reliable genome-wide heteroplasmy map, which can provide further insights into mtDNArelated diseases and the evolution of mtDNA.