Pedigree derived mutation rate across the entire mitochondrial genome of the Norfolk Island population

Estimates of mutation rates for various regions of the human mitochondrial genome (mtGenome) vary widely, depending on whether they are inferred using a phylogenetic approach or obtained directly from pedigrees. Traditionally, only the control region, or small portions of the coding region have been targeted for analysis due to the cost and effort required to produce whole mtGenome Sanger profiles. Here, we report one of the first pedigree derived mutation rates for the entire human mtGenome. The entire mtGenome from 225 individuals originating from Norfolk Island was analysed to estimate the pedigree derived mutation rate and compared against published mutation rates. These individuals were from 45 maternal lineages spanning 345 generational events. Mutation rates for various portions of the mtGenome were calculated. Nine mutations (including two transitions and seven cases of heteroplasmy) were observed, resulting in a rate of 0.058 mutations/site/million years (95% CI 0.031–0.108). These mutation rates are approximately 16 times higher than estimates derived from phylogenetic analysis with heteroplasmy detected in 13 samples (n = 225, 5.8% individuals). Providing one of the first pedigree derived estimates for the entire mtGenome, this study provides a better understanding of human mtGenome evolution and has relevance to many research fields, including medicine, anthropology and forensics.


Results
Whole mtGenome sequencing to identify mtDNA variants in the Norfolk Island pedigree. The whole mtGenome was sequenced for 225 individuals representing 45 independent mtDNA lineages from the core pedigree of the NI population isolate on an Ion Torrent NGS platform (Fig. 1). Sequence quality (Phred) scores remained consistent at > 25 for all samples at the median read length (140 bp), with data reaching a median depth of approximately 370X for all 225 samples. Sequencing reads were aligned in relation to the rCRS and variants annotated using Mitomaster, with called variants subsequently validated by Sanger Sequencing.
In this study, NGS analysis of homopolymer stretches (such as poly C regions) were a challenge for the Ion Torrent platform resulting in sequencing slippage in regions containing long homopolymer stretches. This is due to the degree of change in voltage loss resolution (above 6-8 bp) inherent to the platform. As a result, these regions have poor mapping quality and can lead to false-positive and false-negative calls and typically require assessment using alternative methods. The identification and calling INDELS also presented a bioinformatic challenge in the samples due to sequencing error bias and required a higher level of manual assessment. Issues in sequencing data in homopolymeric regions are not unique to NGS. Sanger sequencing is also susceptible to such errors 24 , with all sequencing platforms prone to sequencing errors corresponding to the chemistry and platforms used 25 .
In the rCRS, position T16189 is flanked by cytosine residues between positions 16183 and 16194. In the samples examined in this study, this sequence stretch was consistently and reproducibly represented in NGS and Sanger data analysis when no length or point variation was present. However, in 134 samples (n = 225, 59.6%), www.nature.com/scientificreports/ a combination of the T16189C transition and a deletion at position 16189 resulted in low coverage values and unreliable Sanger data for this region. Consequently, it became difficult to determine the number of cytosine residues present. In accordance with Interpretation Guidelines for mtDNA Sequencing outlined by the Laboratory Division of the FBI 26 , no attempt was made to count the number of residues for interpretation purposes and all comparisons assumed the same number were present. Similar results were obtained for the HVII C-track between position 302 to 310 and 310 to 316. Difficulties were also observed when an insertion of one or more Cs (resulting in a stretch of ≥ 8 C residues) were present. Mutations in these regions were omitted during the comparison of maternally related individuals, and hence were not included in the mutation rate calculations performed in this study.
Calculation of mutation rates across the mtGenome. When the individual pedigrees were examined, the number of meioses and mutations were observed (Table 1). Homoplasmic mutations were found at two different sites in one individual each. Heteroplasmy was confirmed at 7 positions across 13 individuals (n = 225, 5.78% individuals). Analysis of the individuals belonging to 45 families from NI revealed 13.3% of the families contained at least one individual with mtDNA heteroplasmy as a result of mutation across the entire mtGenome. The frequency of point heteroplasmy in the NI samples did not show significant differences in relation to gender (Chi-squared test: X 2 = 0.472, df = 1, P = 0.492), nor was it found to be associated with a specific mtDNA haplogroup (Fisher exact test: P = 0.080). Further analysis of the 345 mtDNA transmissions identified 9 mutations (2 transitions and 7 cases of heteroplasmy) across the entire mtGenome, suggesting that 1.57 × 10 -6 mutations occur (at a detectable level) in the mtGenome in each generation. Assuming a generation time of 26.9 years, the mutation rate for the entire mtGenome would be 0.058 mutations/site/Myr when including heteroplasmy (95% CI 0.031-0.108), and 0.013 mutations/site/Myr when excluding heteroplasmy (95% CI 0.003-0.047) ( Table 2). In order to make an estimate of separate mutation rates for the HVI and HVII regions, HVI was delimited by positions 16024-16383 (360 bp) and HVII delimited by positions 57-371 (315 bp) in accordance with Sigurgardottir et al. 2 .

Discussion
Mutation rate. In the analysis of 345 genetic transmissions, 9 mutations (2 transitions and 7 cases of heteroplasmy) were detected across the entire mtGenome, suggesting the mutation rate for the entire mtGenome to be 0.058 mutations/site/Myr (95% CI 0.031-0.108). Table 3 compares the derived mutation rate for each region calculated using the NI samples with previous studies using closed populations. NI sample derived mutation rates are shown in bold. For consistency, mutation rates for each study were converted to mutations/site/Myr with a defined generation time of 25 years. Only rates that could be accurately standardised were included in the table. As such, a number of published rates were excluded, for example those by [28][29][30][31][32][33][34][35] .
Note that although rates are grouped based on region (such as control region, coding region), the sequence range varies from study to study. Furthermore, inclusion requirements varied for each study impacting the rate calculated (for example, only substitutions for which there was evidence of a germinal origin were considered in Santos et al. 9 , or all substitutions detected were considered in Madrigal et al. 36 ). As expected, the mutation rate for the coding region was found to be smaller than that reported for the control region, and as with other studies, the estimation reported here is much higher than those obtained by phylogenetic methods (Table 3).
In the NI samples, with the exception of the coding region, the mutation rate calculated for various regions of the mtGenome were found to be relatively high when compared to previous published studies. This may be attributed to differences in the population examined, study design or analysis. In this study, the rates calculated for the NI population included all mutations (heteroplasmy, germline, somatic). In contrast, for the HVI region, Santos et al. 9 reported a mutation rate of 0.074 mutations/site/Myr (16 times lower than the rate for NI: 1.2 mutations/site/Myr-both using a generation time of 25 years), however their rate only included substitutions (including heteroplasmy) with a germinal origin present in women fixed at the individual level. While inclusion criteria alone does not fully explain the variations identified via pedigree-derived mutation rates, it likely contributes to the observed difference.
Interestingly, the mutation rate for the coding region was one of the smallest rates reported, despite our study being the only one to examine the entire region (bases 576-16024). This variation could be due to population differences. NI is a small, isolated population, and therefore the evolution of extremely high mutation rates is unlikely to occur unless organisms are under special circumstances 37 , with beneficial mutations rarely compensating for deleterious mutations.
In addition, heteroplasmic variants were not included in all studies assessing the mutation rate of mtDNA (for example 1 ). It is reasonable to suggest that heteroplasmy may resolve in favour of a 'new' base, and therefore, exclusion of heteroplasmic variants may underestimate the true mutation rate. However, the inclusion of heteroplasmic mutations introduces the limitation that only those that have reached a level of 20% of the mtDNA population are deemed detected. This raises the question of whether the discrepancy between phylogenetic and pedigree derived rates could be answered by simply excluding heteroplasmic variants, on the presumption that they will not become homoplasmic and/or will not be transmitted at the level of the population. Heteroplasmic mutations that reached 20%, the minimum threshold required to be included in the pedigree derived mutation rate reported here, should have a higher probability of becoming homoplasmic than those reaching only 1-2%.

Heteroplasmy.
Several previous studies have demonstrated Sanger sequencing to be valid for quantification of heteroplasmy greater than 10% and that NGS could detect and quantify heteroplasmy as low as 1% 45 . In this study, to ensure robustness of the data and subsequent analysis, a threshold of 20% or greater was used. Previous studies have also performed independent DNA extraction, PCR amplification and sequencing to authen-  www.nature.com/scientificreports/ upon the same method to quantify levels of heteroplasmy. Specifically, these studies measured the height of the two peaks directly from the sequencing electropherogram and calculated the proportion between the heights of each peak with respect to the sum of the height of the two peaks. For comparison purposes, the obtained numerical proportions were then averaged. Another method previously used included cloning PCR products encompassing the D-loop region and subsequent sequencing of the 26-66 clones to determine the number containing each mtDNA variant 9,46 . The increased percentage of heteroplasmy observed in those studies when compared to ours could be due to the lower minor variant threshold used. The threshold used in Ramos et al. 46 was 10%, and although no threshold was stated in Santos et al. 9 , that study included cases of heteroplasmy where the minor variant was as low as 2.5%. In addition to those outlined in Table 1, heteroplasmic mixtures below the 20% threshold were observed in three additional individuals, with MAFs of 12, 15 and 15% for A16280R, T9012Y and C16320, respectively. As these were below the 20% threshold, they were not considered in any further analysis. Further heteroplasmic variants may be evident at other positions across the mtGenome if our threshold was reduced to one comparable with Ramos et al. 46 or Santos et al. 9 . In our study, 6.7% (n = 45) of the NI pedigrees examined, at least one individual within the lineage presented mtDNA heteroplasmy produced by mutations across the coding region. This value is consistent with the percentage reported by Santos et al. 43 (6.5% of families). Interestingly, no individuals from the NI pedigrees examined in this study showed heteroplasmy at the positions examined by Santos et al. 43 .
Length heteroplasmy is typically observed in every individual where a transition at position T16189 results in a homopolymer of nine or more cytosine residues, with no length heteroplasmy observed when seven or fewer cytosine residues are present [48][49][50][51] . In the NI samples, consistent and reliable Sanger sequencing was not achieved for individuals that exhibited the T16189 transition and therefore, length heteroplasmy could not be confirmed. Hence, to ensure consistent reporting, the decision was made to disregard variants in positions 16180-16183 for all samples. The poly-C tracks of both the HVI and HVII regions are known to have high insertion/deletion rates, resulting in length heteroplasmy 52 . For example, Santos et al. 9 reported length heteroplasmy produced by the insertion of cytosine residues in the poly-C tract of HVRI and HVRII respectively in 22.92% and 54.16% of the families studied. It is accepted that the general mechanism for generating length heteroplasmy in these regions is replication slippage 53 . These regions were also excluded during the analysis of the NI individuals and no other length heteroplasmy in any region was observed in the NI samples.

Conclusion
The mtDNA control region is commonly used to assess human evolution and population movements. Many of these studies rely on phylogenetic analysis of mtDNA control region haplotype trees and phylogenetically derived rates of divergence. Estimates of the mutation rate with a non-phylogenetic approach (namely pedigree analysis) are reported to be approximately ten-fold higher than phylogenetically derived rates. Using a pedigree approach, 9 mutations (2 transitions and 7 cases of heteroplasmy) were identified across the entire human mtGenome, resulting in a mutation rate (obtained by employing the same definition of mutation used by other authors) of 0.058 mutations/site/Myr (95% CI 0.031-0.108). The mutation rate produced from this study is one of the first to use extended pedigree analysis of the entire human mtGenome in combination with NGS. If mutations arise randomly in the mtGenome, these results suggest that newly arising mutations in the human mitochondrial coding region are eliminated before they reach a detectable frequency. Defining this mechanism may provide further insight into the interpretation of human mtGenome data across numerous fields of study.

Materials and methods
Ethics. Ethical clearance for the Norfolk Island mitochondrial DNA analysis portion of this study was provided originally by the Griffith University Human Research Ethics Committee (Approval MSC/04/09/HREC). Ethical clearance was transferred to and is now provided by the Queensland University of Technology Human Research Ethics Committee (Approval Number: 1400000749). No other ethical clearance was required. All analyses were performed in accordance with relevant guidelines and regulations, and all Norfolk participants provided informed consent for research involvement.  [17][18][19][20][21] . The NI core pedigree contains individuals that are most closely related to the original founders. For illustration, the NI core pedigree is shown in Fig. 1. Using the NI core pedigree, maternal pedigrees were generated by establishing a list of founding mothers and tracing their maternal lines. Resultant pedigrees were generated using the Pedigree v1.4 54 , Kinship2 v1.8.5 55  Data analysis. Data analysis was performed using the bioinformatics pipeline outlined in Harvey et al. 58 .
Sequences were obtained for the entire mtGenome and were aligned in relation to the revised Cambridge Reference Sequence (rCRS) 59 using an online tool Mitomaster to annotate variants and call haplotype information for each sample 60 .  61 . Variants in positions 16180-16183 were also ignored as reliable sequence data could not be obtained (see "Results" and "Discussion" sections).
In cases of lengthy heteroplasmy, a single dominant variant was identified. In cases of point heteroplasmy, the position was deemed heteroplasmic if the following requirements were achieved: 1. For NGS data, the minor allele frequency was greater than 20% (of total coverage). This threshold was chosen to reliably differentiate between signal and noise and reduce false positives.

Figure 2.
Forty-five maternal pedigrees relating to the sampled individuals examined in this study. Individual pedigrees were established for forty-five founding mothers from the Norfolk Island Core Pedigree (Fig. 1). Individuals whose mtDNA was sequenced are shown as blackened circles (females) or squares (males) and are from the lower four generations of the Norfolk Island Core Pedigree (Fig. 1 www.nature.com/scientificreports/ 2. When sequenced using Sanger sequencing, the minor allele was readily visible (a distinct peak of normal morphology was evident and white space beneath the peak could be observed without changing the chromatogram view to examine the signal closer to the baseline) 48 . 3. The minor allele was evident in two different high-quality Sanger sequences (for example, when using both forward and reverse primers).
All mutations were confirmed using Sanger sequencing. Amplification was performed in a 25 μL total reaction volume using 1.6 μL magnesium chloride, 5 μL 5X GoTaq Colourless Flexi Reaction Buffer, 0.5 μL of 10 mM dNTPs, 9.4 μL deionised water, 2 μL GoTaq Hot Start Polymerase, 2 μL of each 10 μM amplification primer and 50 ng DNA extract. Thermal cycling conditions were 95 °C for 30 s; 30 cycles of 95 °C for 30 s, 57 °C for 40 s, 72 °C for 30 s; and a 68 °C hold for 3 min. Purification of amplified products prior to sequencing was performed using ExoSAP-IT PCR Product Clean-up Reagent, using 5 μL PCR product, 4 μL deionised water, and 1 μL ExoSAP-IT. Thermal cycling conditions were as per the manufacturer's instructions. Sequencing was performed via the BigDye Terminator v3.1 Cycle Sequencing Kit. Reactions consisted of 9.64 μL deionized water; 4 μL BigDye Terminator v3.1 Sequencing Buffer; 0.5 μL BigDye Terminator v3.1 Ready Reaction Mix; 3.2 μL sequencing primer at 10 μM; and 2.66 μL PCR product for a total reaction volume of 20 μL. Thermal cycling conditions used were: 96 °C hold for 1 min; 30 cycles of 96 °C for 10 s, 50 °C for 5 s, and 60 °C for 4 min; followed by 1 cycle of 4 °C for 5 min and 10 °C for 5 min. Sequence product purification was performed via ethanol precipitation. Sequence detection was performed by capillary electrophoresis on a 3500 Genetic Analyser using a 50 cm array, the FastSeq instrument protocol with default instrument settings. Amplification and sequencing primers are outlined in Supplementary Table S1.
The mutation rate was derived from the number of detected mutations per number of 'meioses' or transmission events, which is the number of cumulative generations tracing back to the maternal ancestor. Genealogical records for a sample of 222 individuals from Norfolk Island suggest that the average maternal intergenerational time is 26.9 years. The mutation rate is expressed as mutations per base pair (site) per million years (mutations/ site/Myr), where the generational time was assumed to be 26.9 years. Different values of the mutation rate were estimated according to different assumptions: (1) all mutations that were detected (including heteroplasmy) were considered for the mutation rate calculation, and (2) all mutations that were detected (excluding heteroplasmy) were considered for the mutation rate calculation.
Confidence intervals (CIs: 95%) were calculated using Epitools, an online tool provided by AusVet Animal Health Services 62 . The program outputs intervals using five alternative calculation methods as described in Brown et al. 63 . The Wilson and Clopper Pearson methods were reported.

Data availability
The data that support the findings of this study have restricted access due to ethical requirements and agreements with the Norfolk Island community, and so are not publicly available. Aggregate data may, however, be made available from the authors upon request. The Norfolk genetics steering committee will assess data access requests via our GRC computational genetics group (interested researchers should contact grccomputational-genomics@gmail.com).