A new method for long-read sequencing of animal mitochondrial genomes: application to the identification of equine mitochondrial DNA variants

Abstract

Background

Mitochondrial DNA is remarkably polymorphic. This is why animal geneticists survey mitochondrial genomes variations for fundamental and applied purposes. We present here an approach to sequence whole mitochondrial genomes using nanopore long-read sequencing. Our method relies on the selective elimination of nuclear DNA using an exonuclease treatment and on the amplification of circular mitochondrial DNA using a multiple displacement amplification step.

Results

We optimized each preparative step to obtain a 100 million-fold enrichment of horse mitochondrial DNA relative to nuclear DNA. We sequenced these amplified mitochondrial DNA using nanopore sequencing technology and obtained mitochondrial DNA reads that represented up to half of the sequencing output. The sequence reads were 2.3 kb of mean length and provided an even coverage of the mitochondrial genome. Long-reads spanning half or more of the whole mtDNA provided a coverage that varied between 118X and 488X. We evaluated SNPs identified using these long-reads by Sanger sequencing as ground truth and found a precision of 100.0%; a recall of 93.1% and a F1-score of 0.964 using the Twilight horse mtDNA reference. The choice of the mtDNA reference impacted variant calling efficiency with F1-scores varying between 0.947 and 0.964.

Conclusions

Our method to amplify mtDNA and to sequence it using the nanopore technology is usable for mitochondrial DNA variant analysis. With minor modifications, this approach could easily be applied to other large circular DNA molecules.

Background

Mitochondria have become crucial organelles of eukaryote cells since an ancestral bacterial endosymbiosis event gave rise to the eukaryote branch of the tree of life, about 1.5 billion years ago [1, 2]. The discovery that mitochondria contains DNA in the early 1960s attracted the attention of many scientists and more recently the word of mitogenomics has been coined to unite this field of study [3]. Today, thousands of mitochondrial genomes (mtDNA) from many species have been sequenced, and ongoing studies require robust, rapid and cheap method to decipher novel genomes and identify genetic variants [4]. In human alone, more than 50 thousand full length mtDNA sequences are known, but much less sequences and variants are available for other mammals [5]. Numerous and diverse applications take advantage of mtDNA sequencing, such as diagnostic of mitochondrial diseases in human and animals, breeding selection, forensic studies, biodiversity surveys, conservation genetics and evolutionary analysis.

Animal mitochondrial genomes are usually encoded by a single circular DNA of about 10 to 20 kbp, but weird cases are known such as linear or very large molecules [6]. The mitochondrial genomes are characterized by a maternal transmission, a lack of recombination and higher frequency of mutations than in the nuclear genome [7]. In vertebrates, the ratio of mtDNA over nuclear DNA mutation rate is usually above 20 [8]. This means that mtDNA polymorphisms may occur relatively often in animal populations with the potential for associated disorders, and this is why sequencing these organelle genomes is an important component of contemporaneous animal genetics. Knowing the whole mtDNA sequence enables the identification of all possible genetic variants: single-nucleotide polymorphisms (SNP), insertions, deletions, duplications, inversions, rearrangements. Although mitochondrial variants are typically maternally inherited, they can be sporadic. Such de novo variants are difficult to assess due to the co-existence of several mtDNA population in a cell (heteroplasmy) that may be below an assay’s detection level, or different between tissues of the same organism [9]. The presence of DNA segments transferred from the mitochondria to the nuclear DNA, the so-called Nuclear Mitochondrial DNA Sequences (NUMTs), can seriously complicate mtDNA analysis in a variety of animals from bees to humans [10,11,12] . Such NUMTs can span several kilobases, can be nearly identical to mtDNA sequence and can exhibit variation between individuals of the same species [13].

The current methods for mtDNA sequencing rely on standard sequencing methods such as Sanger sequencing or massively parallel sequencing [4]. The choice of the starting DNA material is of importance since the relative mass amount of mtDNA is very low in comparison to nuclear DNA. Typical DNA preparations comprise approximately 0.1% mtDNA, even though there are hundreds to tens of thousands copies of mtDNA per cell [14]. One strategy relies on the amplification of mtDNA using regular or long-range PCR. This process is associated with an unpredictable risk of PCR errors, it can introduce a sequence bias, and it can be tricky because it depends on primers, on PCR efficiency, and finally the whole procedure can be time consuming and costly. In addition, the reliance on PCR primers derived from a reference sequence can drive the failure of some amplifications because they may suffer from a lack of specificity to mtDNA, and they may lead to the production of artefactual amplicons, including some derived from NUMTs [15,16,17]. The other strategy rely on whole-genome shotgun using massively parallel sequencing approaches, an efficient approach for mitogenomics because mtDNA reads are more prevalent than nuclear DNA reads, but it is not bullet-proof and short read mapping can also be confounded by NUMTs [18, 19]. These approaches involve significant resources for bench work and efforts in sequence handling, they are cost-effective only when high number of samples are pooled and therefore they do not scale down easily [18]. The principal caveat of the current mtDNA sequencing methods is that raw sequences enable to read only a tiny portion of the DNA at a time, preventing definitive evidence for long-range phasing or structural variant identification.

Long-read sequencing using Oxford Nanopore or PacBio technologies provide the possibility to capture in a single read the majority or the whole sequence of mtDNA molecules. Therefore, it can provide definitive evidence that a given mitogenome exists in a single cell or in a sample. So far, a few studies harnessed the Nanopore technology for sequencing mtDNA. Istace and colleagues assembled 21 yeast mitogenomes in a study where they evaluated the performances of nanopore-only assemblers. They used only 2D reads for mitogenome assembly, a chemistry that provided better sequence quality but that has been abandoned [20]. Ranjard and collaborators assembled the green-lipped mussel mitogenome (Perna canaliculus Gmelin 1791) using a whole-genome shotgun approach with long-reads to obtain a draft assembly that was polished using short-reads from a RNAseq experiment; but they do not give any details on the long-read sequencing output [21]. In another study, the european lobster (Homarus gammarus) mitogenome was assembled separately with both Illumina short-reads and Nanopore long-reads [22]. The report by Gan and collaborators is instructive since they could not obtain a large mtDNA contig using de novo or single gene bait-based assemblies of six million Illumina paired-end reads totalling 1.5 Gbp. However, the 20 kbp lobster mitogenome could be assembled from 40,687 long-reads totalling 100 Mbp. The same authors reported another mitogenome assembled similarly, with only nine long-reads [23]. These reports testify the usefulness of the Nanopore technology to sequence mitochondrial genomes, but the corresponding protocols are not tailored to identify mtDNA variants on a large-scale.

In horses, many mitochondrial functions have been shown to play important roles for muscular exercise and disorders [24,25,26,27]. Yet, the study of mtDNA variation is lagging behind even though it has been instrumental to study ancient horse’s DNA and unravel the domestication events of the wild horse [28,29,30]. As in other mammals, the abundance and polymorphism of NUMTs is a characteristic of the horse genome and this has negative impact on the reliability of mtDNA variant identification [12, 31, 32].

The aim of our study was to trial a new method for mtDNA resequencing using the long-read technology provided by Oxford Nanopore Technologies (ONT) with the goal of identifying mtDNA variants. We developed an optimised method for mtDNA sequencing starting from total genomic DNA obtained using routine extraction methods. We implemented a protocol based on a nuclear DNA removal by nuclease treatment followed by mtDNA enrichment using multiple displacement amplification (MDA) and sequencing using the 1D chemistry on the MinION platform of ONT. We anticipate that this method will ease the identification of mtDNA variation.

Results

Enrichment of circular mtDNA by exonuclease V treatment

We started our strategy by enriching the circular forms of animal mtDNA in genomic DNA extracts. Since traditional biochemical purification of mitochondria does not scale up easily and requires a large amount of fresh tissue, we selected an alternative method. We chose an exonuclease treatment to deplete linear DNA fragments and thereby improve the proportion of mtDNA versus nuclear DNA, as previously described [33]. A frequently used enzyme to degrade linear DNA molecules is the exonuclease V that has different nuclease activities, comprising an ATP-dependent double-stranded and bi-directional exonuclease activity. Since the amount of nuDNA vastly exceeds that of mtDNA, we tried four conditions of nuclear DNA removal and evaluated their efficiency using qPCR (Fig. 1a). The effect of exonuclease V treatments impacted dramatically the quantity of nuDNA, since it diminished by several thousand fold relative to mtDNA. mtDNA quantity diminished slightly after the various exonuclease treatment to only about half of its amount in the starting DNA extract. We observed that increasing the incubation time had an effect while the amount of exonuclease V seemed not to be a limiting factor. We found that the best condition was 2 h of incubation with 10 Units of Exonuclease V, since it resulted in a ratio of mtDNA over nuDNA of 3755+/− 0.7 (Fig. 1a).

Fig. 1

Quantification of mtDNA enrichment by qPCR. a Quantification after nuclear DNA removal using exonuclease V (ExoV) as indicated. b Quantification after mtDNA amplification using MDA. The results are given for DNA untreated by exoV (No ExoV), or after two conditions of ExoV digestion as indicated. REPLI-g was performed using primers at a concentration of 0.2 μM or 1.0 μM

Such an enrichment was of interest, but it came at the cost of a limited amount of DNA which is inadequate for sequencing without an amplification step. We decided to amplify the mtDNA using a Whole Mitochondrial Genome Amplification (WMGA) method based on multiple displacement amplification [34]. We optimized primer design and concentration in the MDA reaction to avoid template-independent amplification, a technical bias usually observed with MDA. To evaluate the impact of the Exonuclease V treatment and the amplification efficiency of WMGA, we estimated the relative amounts of mtDNA versus nuDNA with or without an exonuclease V treatment, and using two primer concentrations. We present our results on the quantification of the amplification efficiency in Fig. 1b. We observed a nearly 1000 fold better amplification using WMGA on DNA treated with exonuclease V versus the starting undigested DNA (Fig. 1b), with an overall amplification of more than 100 million fold in the best condition using 0.2 μM of each primer. We obtained on the order of ten micrograms of DNA after WMGA, with an average fragment size of about 3 kbp that we deemed usable for long-read sequencing (SI_Figure 1).

Long-read sequencing

We performed a first sequencing run on nine DNA samples obtained from three horses (H25, H26 and H27) after WMGA on DNA treated or not by exonuclease V and amplified using the two primer concentration conditions mentioned previously. We will refer to DNA samples untreated by exonuclease and amplified using the lowest primer concentration as “exo –“; DNA samples treated by exonuclease and amplified using the lowest primer concentration as “exoA” and DNA samples treated by exonuclease and amplified using the highest primer concentration as “exoB” (Fig. 2, Table 1). We obtained a total of 634,062 reads with a mean length of 2370 nt and a mean quality of 10.2. After demultiplexing, we obtained 562,981 reads (88.8%) evenly distributed among the nine barcodes; with between 45,641 reads and up to 67,448 reads (Table 1 and SI Table 1). We did not observe gross differences of read size or quality between the nine barcoded libraries (Fig. 2, Table 1). We mapped these reads to the reference horse genome and counted the reads mapping to the mitochondrial genome and those mapping to the nuclear genome (Fig. 2b). The lowest proportion of mtDNA reads, about 1%, was from the exoB DNA samples (Fig. 2b). We observed around 5% of mtDNA reads in the exo- DNA samples (Fig. 2). The proportion of mtDNA reads was up to ten times higher in the exoA DNA samples (Fig. 2), with values ranging from 17,569 to 40,888 reads, i.e. 30 to 59% of the total sequences from a given sample. Upon inspection of the location of reads mapping to the nuclear genome, we found that about 20% of them were mapped to NUMTs. The effect of the exonuclease treatment and the amplification on the mapping of reads to mtDNA versus nuDNA was significant (Chi2 = 12,552.68, p < 0.001). The distribution of read size and quality was similar whether they mapped to the mitochondrial or to the nuclear genome (SI Figure 2). In conclusion, the protocol we developed to amplify mtDNA led to a significant increase in the proportion of long sequences derived from the horse mitochondrial genome.

Fig. 2

Read length distribution and proportions of sequences mapping to the mitochondrial or nuclear genomes. a Distribution of reads length mapping to mtDNA (top) or nuDNA (bottom). For clarity, values are shown for a randomly selected subset representing 0.005% of reads. b Numbers and proportions of sequences mapping to mtDNA or nuDNA. The number of mapped sequences is shown on the y axis. exo-: no Exonuclease V treatment; exoA: Exonuclease V 10 U 2 h + REPLI-g 0.2 μM primers; exoB: Exonuclease V 10 U 2 h + REPLI-g 1.0 μM primers

Table 1 Sequencing statistics. H25, H26 and H26 are the three different horse’s DNA analysed. Barcode01 to Barcode09 are the names of the barcodes used in the multiplexing with the kit EXP-NBD-104. The category Unclassified refers to reads that were filtered out at the demultiplexing steps because the barcode sequence was not attributed. exo -: no Exonuclease V treatment; exo A: Exonuclease V 10 U 2 h + REPLI-g 0.2 μM primers; exo B: Exonuclease V 10 U 2 h + REPLI-g 1.0 μM primers

Overall coverage of the mitochondrial genome

We then looked at the coverage of the mitochondrial genome at the base level (Fig. 3). Whatever the treatment, we found that the coverage of the mitochondrial genome was complete and this pattern was reproducible between the samples. We could observe the effect of some primers on the coverage plot that appeared as bulges, for example between positions 13,946 and 14,748. For samples untreated by exonuclease in the S1 sequencing run, the minimum coverage was 66 reads at a given mitochondrial base (SI Table 1). In the same run and for the samples treated with exoA the minimum coverage was 526 reads. The variation between minimal and maximal coverage was between 2.67 (sample H26 exo-) and 5.33 (sample H27 exoB). In a second sequencing run, the coverage was more homogeneous between DNA samples with variations of 3.26–3.97 between minimal and maximal coverage, a minimum coverage of 936 and a maximum of 11,707 (SI Figure 3). In conclusion, the read coverage of the mitochondrial genome was complete and unbiased in our tested conditions.

Fig. 3

Sequence coverage of the mitochondrial genome. The base position coverage is shown on the y axis. Sequences derived from the exo minus samples are shown in light grey, exo plus A are in dark grey and exo plus B are in black. Coverage was computed from mapping sequences to a duplicated mitochondrial genome. The grey vertical line at position 16,666 represents the end of the reference mtDNA sequence. Dotted vertical lines indicate the starting positions of the primers used for MDA

Long-read distribution and coverage of the mitochondrial genome

Since we obtained a good coverage of the mitochondrial genome, we wondered about the proportion of long-reads covering a significant portion of the mitochondrial DNA. In H25, H26 and H27 exoA samples, we found respectively 52, 144 and 116 reads longer than half the mitochondrial genome while we only found between one and seven such reads in the exo- and exoB samples (Fig. 4, SI table 1). In about 10% of these long-reads, we observed that the longest alignment between the read and the mitochondrial genome was remarkably smaller than the read length even though we selected reads for which the alignment coverage on mtDNA exceeded 80% (Fig. 4a). Upon inspection of these cases, we identified two categories of such reads. The first category corresponded to reads spanning more than one full revolution of the circular mtDNA. For example, we obtained five reads longer than the mitochondrial genome, one of which measuring 17,501 nt and spanning more than one full turn with about 1 kbp of additional coverage (Fig. 4). The second category of reads were usually compound and made of two inverted repeated sequences. These could be reads originating from the sequencing of both strands of the same molecule, one after the other like in the 2D or 1D2 approaches developed by ONT, with usually the second strand being incomplete. Such reads could correspond to non-chimeric or chimeric DNA molecules produced during WMGA. Altogether, these long-reads covered the three horse’s mitochondrial genome at 30.8X, 87.8X and 70.1X (Fig. 4b). In conclusion, we obtained very long-reads from the exonuclease treated DNA spanning large portions of the mitochondrial genome up to its entirety.

Fig. 4

Mitochondrial long-reads. a Scatter plot of alignment and read lengths. Dotted lines represent the size of the mitochondrial genome. Each square, triangle or dot represents one read satisfying the following criteria: the read aligns to more than 8000 bp of the mtDNA reference sequence and more than 80% of the read aligns to the reference mtDNA sequence. Circles, triangles and squares label different read categories correspond to the ratio of read length over alignment length: « normal » for values less than 1.1, « middle » for values between 1.1 and 1.2, « high » for values over 1.2. Values higher than 1 mean that there are mutiple alignments between the read and the mitochondrial genome; possibly because the read spans more than one full turn of the mtDNA when the read is longer than the mtDNA or when it is « 2D like ». Most of these reads align as expected over their whole length to a single portion of the mtDNA reference sequence (circles that can be seen along the diagonal). A small portion of these reads (triangles and squares) show more complicated alignment figures with the mtDNA reference sequence. We called some reads « 2D like » because their sequence encompass the same mtDNA segment on the same or opposite orientation. b Long-read coverage display on the horse mitochondrial genome. Each horizontal grey bar corresponds to one read satisfying the criteria defined in a, and is displayed according to its length on the y-axis. The read matching mtDNA region is indicated by the span of the bar on the x-axis (mtDNA bp coordinates). Remark the longest read (17.5 kbp) that spans more than the length of the mtDNA with the region covered twice highlighted by the dark grey color

These results of the first sequencing run prompted us to reproduce in a second set of experiments the analysis of the same set of horse DNAs. We prepared sequencing libraries using the LSK109 chemistry from 12 new samples treated using the exoA conditions. We changed the primer pool composition, excluding three primers matching to numerous NUMTs loci in the nuclear genome, performed the exonuclease treatment in duplicate and made a technical replicate of each experiment to have a better readout of the factors that could affect the efficiency and the reproductibility of our protocol. We multiplexed the 12 resulting samples for sequencing on the same flow cell. We obtained more reads from this second sequencing run and therefore a much higher coverage of the mitochondrial genome (SI table 1). The fraction of mitochondrial reads ranged from 16 to 20%, with the exception of two technical replicates reaching 50%. Like previously, the mitochondrial coverage was unbiased (SI Figure 3). The mitochondrial genome coverage by long-reads (> 8000 nt) varied between 118X and 488X, and because the sequencing yield by sample was higher, we also obtained significant coverage by reads longer than 16 kb, from 2 to 40X (SI Table 1, SI Figure 4). In conclusion, we obtained reproducible results in terms of mitochondrial reads enrichment, mitochondrial genome coverage and long-reads.

Mitochondrial DNA variants calling

We then performed variants calling on all sequencing data sets obtained from the three different horses DNA, i.e. from the seven sequencing libraries (three from the first sequencing experiment and four from the second sequencing experiment). For this purpose, we used sequence reads of 2000 nt or longer and satisfying a minimal quality threshold of 10, since we obtained a deep coverage and we wanted to avoid reads potentially derived from NUMTs. We used the medaka variant calling method and identified a total of 88 variable positions based on the same Arabian breed mtDNA reference as before: 77 for H25 samples, 10 for H26 and 83 for H27 (Fig. 5a). All variants were transitions, and all were found in at least one known horse mtDNA sequence available in GenBank. The majority of these variable positions were identified in the seven sequencing libraries prepared from the same DNA: 73/77 for H25 (94.8%), 10/10 for H26 (100.0%) and 77/83 for H27 (92.8%). There were 52 coding and 36 noncoding SNPs, 26 of which from the control region. The remaining positions were usually of low quality and were counted as false negatives, because they were missed in at least one sequence sample (see below).

Fig. 5

Variant calling from long mtDNA reads. a The top panel represents the mitochondrial genome with the scale in base pairs and its rRNA, tRNA, gene and control region features. The three mtDNA reference sequences were aligned to visualise their differences. Remark the significant differences in the control region. The position of the identified polymorphisms in the H25, H26 and H27 exoA samples are indicated just below as small black vertical lines. b The lower panel is an heatmap showing all variants identified along the mtDNA (from the left to the right) for the different samples indicated and the three different references (from top to bottom). Variants are coloured according to their quality score or are shown in white when they were not found in a given horse DNA. SNPs categories defined by comparison with Sanger sequencing are shown below the heatmap: false negatives (FN, in red), false positives (FP, in orange) and true positives (TP, in green). Coding SNPs are shown in black in the last row, while non-coding SNPs are shown in grey

We reasoned that a reference allele bias effect could impact the quality of variant calling, especially in the highly polymorphic and repetitive control region [35]. Therefore, we performed the same variant calling procedure using two alternative references: a Thoroughbred breed more closely related to H25 and H27, and the Twilight horse (a Thoroughbred) used for the equine reference genome. The results are presented in Fig. 5 and details are provided in SI Table 2. We identified 82 variants when we used the Thoroughbred reference mtDNA sequence: 4 for H25, 74 for H26 and 9 for H27 (Fig. 5b). Similarly, most variants were found in all seven sequencing libraries from the same horse: 2/4 for H25 (50.0%), 72/74 for H26 (97.3%) and 7/9 for H27 (77.8%). We found 50 coding and 32 noncoding SNPs, including 23 in the control region (Fig. 5b). We identified 82 variants when we used the Twilight reference mtDNA sequence: 66 for H25, 85 for H26 and 70 for H27 (Fig. 5b). Again, most variants were found in all seven sequencing libraries from the same horse: 63/66 for H25 (95.4%), 82/85 for H26 (96.5%) and 66/70 for H27 (94.3%). We found 50 coding and 44 noncoding SNPs, including 31 in the control region (Fig. 5b). As expected, the distribution of variants was globally the opposite between the one found using an Arabian breed mtDNA reference and the one found using the Thoroughbred, but the latter was characterized by higher variant qualities (compare qualities on Fig. 5b). The distribution of variants found using Twilight as a reference was more even between the three horses. We found the same relative abundances of polymorphisms in coding and non-coding sequences, and as expected the density of polymorphisms in the control region was the highest.

We evaluated the validity of these polymorphisms by sequencing nine PCR fragments using Sanger sequencing as ground truth. We could thus compare 7493 bp read using Sanger and Nanopore sequencing. The overall sequence identity between the Sanger and the Nanopore sequences was 99.89%. We made a concordance analysis based on variants encompassed by both Sanger and Nanopore sequencing (Table 2 and SI Table 2). The values were affected mainly by the number of variants and by the reference being used, and the results obtained using Twilight are globally more meaningful. Overall, we can underscore that the precision was excellent with a single case of false positive observed on H27 with the Arabian and the Thoroughbred reference. We looked at each sample individually to further analyse the between sample variation of true positive and false negative variant calling (SI Table 2). Overall, the mean true positive rate were 93.42, 93.33 and 94.75% with a false negative rate of 5.71, 5.00 and 5.02% using Arabian, Thoroughbred and Twilight references, respectively. These values underscored the critical impact of the reference sequence on variant calling. We observed that all cases classified as false negatives because they were missed in at least one sample were detected in at least another sample from the same horse DNA and according to the same quality criteria (Fig. 5). There were only two variations that were always missed using Nanopore reads. The first, detected using the Arabian reference, was a short deletion in the control region at position 16,402 in an array of cytosines. The second, detected using the Thoroughbred reference, was a dinucleotide polymorphism in a minisatellite repeat of the control region at position 16,125-16,126. This meant that there is no fundamental technical limitation to detect these variants using currently available nanopore long read sequencing, it is just a matter of thresholds and reproducibility in sequence quality in these regions.

Table 2 Concordance analysis

We studied in more details the polymorphisms of the control region. While seven variants turned as false negatives or false positives using the Arabian reference, there was only one difficult position identified using the Thoroughbred reference (Fig. 5). This difficulty lies at the 5′ end of a minisatellite segment made of 29 GTGCACCT repeats. Upon close inspection, both Sanger and Nanopore sequencing confirmed that the first motif of the minisatellite varied (GTGCACCT- > GCACACCT). Yet there were many long reads that aligned with deletions in the minisatellite region in a reference-dependent manner, leading to lowest prevalence of the correct sequence and to an absence or a low quality variant call (SI Figure 5).

In conclusion, the concordance for SNP calling was very good and depended mostly on the reference used and on variations in the mtDNA control region, a notoriously difficult region to sequence.

Discussion

We combined two previously described approaches to enrich the fraction of sequence reads derived from mtDNA, namely nuclear DNA removal using an exonuclease treatment and amplification using MDA [33, 34]. We found that these two techniques require optimisation steps to obtain a significant mtDNA selective amplification factor. For nuclear DNA removal, the exonuclease V amount and incubation time are likely to be dependent on the source and the quality of the input DNA. Using this step, we obtained a 3000 fold enrichment of the mtDNA/nuDNA ratio, so that theoretically the relative mass amount of mtDNA in the whole cellular DNA increased from 0.63% (considering 1000 mitochondrial genome per cell) to 99.99%. This value is in line with the report of Jayaprakash and colleagues [33].

We obtained on the order of a 100,000 fold enrichment after MDA on DNA untreated by exonuclease V, i.e. about ten times more than the value of 100 to 10,000 reported by Marquis and colleagues [34]. We used the same reagents as Marquis and colleagues, with the exceptions of primers and DNA, and the two primer concentrations we used in REPLI-g lead to similar enrichment. This difference of a factor of ten may be explained by a difference between human and horse mtDNA/nuDNA ratio, by the qPCR quantification method and by the DNA isolation procedure, the latter being known to impact the quantification of mtDNA copy number [36]. We observed that there were less mtDNA reads in samples treated using the exoB condition (1 μM of primers) versus the exoA condition (0.2 μM of primers). A possible explanation is that higher primer concentration could lead to non-specific amplification, including template independent amplification, a well-known issue in MDA [37].

A concern related to MDA is its propensity to produce chimeric molecules, an artefact found in 3 to 6% of reads derived from Illumina sequencing [38, 39]. We think that it is unlikely that chimeric reads affect variant calling since chimerism frequency is more than ten times lower than the prevalence of variants. Moreover, such chimeric reads were found to be mostly formed during Illumina sequencing library constructions [38].

We obtained a ratio of mtDNA to nuDNA reads of up to 59%, but this varied for the same starting DNA between different MDA. Still, we obtained a minimum of 16% of reads derived from mtDNA, and this provided an even coverage of the mitochondrial genome. These results are similar to what is reported using some commercial solutions but less than those obtained using PCR on mtDNA extracted using a plasmid preparation kit [40]. The principal difference here lies in the horse NUMTs, because most reads mapped to nuDNA originated from NUMTs.

We used a set of primer to amplify horse mtDNA that were uniformly distributed along the mitochondrial genome, and we did not observe a significant effect on the coverage. The absence of a strong coverage bias lead by mtDNA amplification using MDA was also found by Marquis and collaborators [34]. In the future, it could be desirable to reduce the number of primers because this could enable the production of larger DNA molecules obtained from MDA. Alternatively, random primers could be used for de novo mtDNA sequencing like in the work of Simison and colleagues [41].

We obtained a sequencing depth of at least 1000X that should allow a much higher level of multiplexing. At the time of initial writing, multiplexing by ligation was limited to 24 but the possibility to multiplex 96 samples using the same protocol is now a reality. According to Quail and colleagues, heterozygous variant calling with a 15x coverage using Ion Torrent or various Illumina datasets was comparable to a coverage of 190x for PacBio datasets [42]. We can deduce conservatively that a similar coverage using Nanopore datasets would be enough to call mitochondrial variants in an efficient manner. This means that we could multiplex four or five times more samples, i.e. a hundred, in a single run. With such a capacity, this method could become very competitive in terms of rapidity, scalability and cost-effectiveness.

We estimated at 99.9% the precision of the obtained sequence, which is satisfactory given that the read median identity was estimated at ~ 89% for 1D sequencing using R9.0 chemistry. The process of MDA is known to be highly accurate and can not be considered as the limiting factor when it comes to long-read sequencing. But at the same time we observed a variant false negative rate of ~ 5% that stems from nanopore sequencing imprecisions around homopolymers of small sizes or around micro or minisatellite segments. The precision is likely to increase through improvement of the sequencing chemistry, like the recent introduction of the R10 flow-cells that combines different pores, and by the improvement of the algorithms used for base-calling. We detected only SNPs using our workflow, and the current level of sequencing precision on R9.4 flowcells is likely to preclude the identification of small indels. Further development of variant calling pipelines may alleviate this limitation, and specific tools to detect structural variants could leverage on the added value of long-reads. Another limitation is that the medaka pipeline is tailored for diploid genomes, but suitable models for variant calling on haploid such as FreeBayes could be used.

The overall quality of variant calling was improved when we used an alternative reference sequence to the point where we did not observe false negatives or false positives in coding regions. All the variants we found were transitions (A- > G or C- > T) and it is known that transitions are dominant among single nucleotide polymorphisms. This effect is more pronounced in mitochondrial genomes but also prevalent in nuclear genomes, both in coding and non-coding regions [43, 44]. The initial explanations advanced to explain this observation were that transitions are more conservative changes, but it seems to be a matter of discussion [45]. Nevertheless, the transition to transversion ratio is a useful tool for quality control in variant calling [46].

It turned out that the single variants that could be missed were in the control region, in a particularly difficult region containing a minisatellite. We found that the main factor involving variant calling accuracy was the choice of the reference used for mapping the reads. We can foresee that the optimisation of mapping parameters in conjunction with the filtering of alignments could dramatically improve the robustness of variant calling [47, 48]. Actually, some level of heteroplasmy involving this minisatellite sequence like previously described in rabbit and other species could explain the lower level of quality of variant sequences in this region [49].

Previous studies have reported the use of Nanopore sequencing to sequence mitochondrial genomes, usually by sequencing long amplicons [50, 51]. In a recent study, Lindberg and colleagues evaluated variant allele frequencies from single or mixed DNA sources using both Illumina and Nanopore sequencing of two overlapping amplicons. In comparison to our findings, they reported a higher precision and recall for single DNA samples [52]. Similarly, the added value of long-read sequencing of amplicons for sorting out complex heteroplasmic conditions was leveraged by White and colleagues [50]. Only few reports mention in some details the prevalence of mtDNA reads in whole genome shotgun nanopore sequencing. In their report on Nanopore DNA sequencing aboard the International Space Station, Castro-Wallace and colleagues found that 29 out of 83,332 reads (0.03%) from mouse DNA were derived from the mouse mtDNA [53]. In another study, Zascavage and colleagues relied on Nanopore long-read sequencing for mtDNA analysis [54]. They reported an average coverage of mtDNA reads that varied between 15 and 118 fold from a whole genome sequencing analysis. In their study, the lower coverage limit was inadequate for accurate variant calling, and the coverage range was too large to warrant a routine use.

Conclusions

Overall, we provided a detailed and benchmarked method that should allow mtDNA genome variant identifications not only for horse geneticists but for other animal geneticists.

The identification of causal mtDNA polymorphisms require accurate variant discovery methods. The capacity to phase mtDNA variants, to identify structural variations and nucleotide modifications is a major asset of Nanopore long-read sequencing. The possibility to do so rapidly and at a relatively modest cost is of added interest, and may be of interest for field studies since the MinION sequencer is portable and able to perform real-time sequencing and base-calling [55]. At present, the technology is still hindered by the limited level of sequence accuracy but this is rapidly changing as sequencing chemistry and base-calling algorithms are evolving [56, 57]. The variety of flow cells and equipment available enables to better tailor the sequencing tools according to the number of samples. We believe that this methodology is easily amenable to sequence circular DNA molecules of similar sizes and even beyond, such as plasts, viral circular DNA or endogenous circular DNA [58].

Methods

Nuclear DNA depletion

We used DNA samples obtained in the course of a previous project and kindly provided by Dr. Céline Robert [59]. The DNA were initially extracted from blood samples of H25 (female, Selle Français), H26 (female, Thoroughbred) and H27 (male, Selle Français). We tested nuclear DNA removal by starting with 1 μg of total DNA using either 10 or 20 Units of Exonuclease V (New England Biolabs) and an incubation of 30 min or 2 h at 37 °C. We performed the enzymatic digestions in 30 μL reactions containing 1 mM ATP in NEB Buffer 4. We stopped the reaction by adding 0.8 μl of 0.5 M EDTA followed by heat inactivation at 70 °C for 30 min. DNA was recovered by precipitation with 4 μl of 3 M NaCl, 2 μl of 5 mg/ml glycogen and 150 μl of 100% ethanol. DNA was finally resuspended in 40 μl of nuclease-free water and stored at − 20 °C until used.

Whole mitochondrial genome amplification (WMGA)

We used the REPLI-g mitochondrial DNA kit (Qiagen) to perform WMGA according to the manufacturer’s instruction with the following modifications. We designed 18 primers using CLC Genomic Workbench v6 (CLC bio) based on a horse mitochondrial genome (GenBank KX669268.1; see Table 1). We adopted the following recommendations for primer design: a length of 10–14 bases, the same number of primers hybridize to each strand; primer hybridization sites are evenly spaced and primers include phosphorothioate modifications between the last three bases of their 3′ end. All primers were purchased from Eurofins Genomics. We tested two primer concentrations: 1 μM and 0.2 μM of each oligonucleotide. The reaction set-up was as follows: 4 μl of DNA treated with Exonuclease V was loaded into a 0.2 ml PCR tube containing 16 μl of RNase DNase free water. 27 μl of REPLI-g mt reaction buffer and 10 μl of horse mitochondrial primer mix (1 μM or 0.2 μM each primer) was added to the DNA and the mixture was incubated at 75 °C for 5 min. We added 1 μl of REPLI-g Midi polymerase to the DNA and incubated at 33 °C for 8 h. Inactivation was done at 65 °C for 3 min. The concentration of amplified DNA was quantitated with Quant-iT™ PicoGreen dsDNA reagent (Molecular Probes) on an Infinite M200 Pro (Tecan). We followed published recommendations to avoid possible contaminations [60].

qPCR experiments

We used qPCR to estimate the ratio of mtDNA over nuclear DNA (nuDNA). The amount of mtDNA was evaluated using three amplicons from the Cox1, Cox3 and Nd2 genes. We also used three nuclear gene loci (Epo, Myod1 and Gapdh) to quantify nuDNA amounts. All primers are given in SI Table 3. As recommended in the REPLI-g mitochondrial DNA kit handbook, the amplified DNA was diluted 1:1000 prior to qPCR. We made two technical replicates for each biological sample, and the same qPCR run contained all samples for a given gene to minimize inter-experimental variations. The 20 μL qPCR reactions were assembled manually. Reactions included 5 μL diluted DNA, 10 μL 2X powered Sybrgreen PCR master mix (Applied Biosystems) and 5 μL primers mix (300 nM each primer). Reactions were incubated in a 96 well optical plate (Applied Biosystems) at 95 °C for 10 min, followed by 40 cycles at 95 °C for 15 s and 60 °C for 1 min in a QuantStudio 12 K Flex instrument (Applied Biosystems).

The qPCR data were analysed using the QuantStudio 12 K Flex software (Applied Biosystems). Cq values were obtained using the auto baseline option and the threshold was manually set to the same value of ΔRn = 0.1 for all amplicons. The Cq values were considered valid if the duplicate measures differed by less than 0.5 cycles, undetermined values were set equivalent to a Ct value of 40. We computed the geometric mean of the Cq values of the three mtDNA and nuDNA amplicons, and then we computed the corresponding fold changes. Formation of primer-dimer and amplification specificity were assessed by melting curve analysis.

Long-read sequencing library preparation

Prior to the library construction, we purified the input DNA samples using a 0.5 volume ratio of Agencourt® AMPure® XP magnetic beads. We quantified DNA using a Qubit dsDNA BR assay and performed a debranching treatment on 8–10 μg of DNA in 100 μL reactions containing 50 Units of endonuclease T7E1 (New England Biolabs) in 1X NEB buffer 2, with an incubation of 2 h at 37 °C. We purified the debranched DNA using one volume of Agencourt® AMPure® XP magnetic beads and checked DNA concentration using a Qubit dsDNA BR assay and DNA integrity of using electrophoresis on a 12 K DNA chip (Experion, Bio-Rad).

These DNA samples were then prepared for sequencing according to native barcoding expansion (kit EXP-NBD-103, ONT) and the 1D Native barcoding genomic DNA protocol (kit SQK-LSK108 for the S1 run and SQK-LSK109 for the S2 run, ONT) with some modifications. We treated 2 μg debranched DNA using the UltraII end repair/dA-tailing enzyme mix (New England Biolabs) followed by a purification using a 0.6 volume ratio of AMPure beads. Native barcodes were added by ligation using the Blunt/TA ligase master mix (New England Biolabs) and the ligation reaction was purified using 1 volume ratio of AMPure beads. Equimolar amounts of each barcoded sample were pooled before adapter ligation using the NEBNext Quick ligation Module (New England Biolabs) and a final purification using a 0.4 volume ratio of AMPure beads with the S Fragment buffer.

For the S1 sequencing run, we used DNA amplified with the 18 primers and multiplexed nine samples using barcodes one to nine. We loaded 195 ng of sequencing library. In this S1 run, each of the three horse’s DNA libraries was derived from three different conditions of exonuclease and amplification, so that each horse’s DNA was independently sequenced three times. For the S2 sequencing run, we used DNA amplified with 15 primers, excluding three that were found to match NUMTs loci, and multiplexed 12 samples using barcodes (barcodes one to twelve). We loaded 439 ng of the sequencing library. In this S2 run, we performed each exonuclease treatment in duplicate and each amplification in duplicate, so that each horse’s DNA was independently sequenced four times. Altogether, each horse’s DNA was sequenced across seven libraries.

We performed the sequencing using R9.4 flowcells from two different batches on a MinION sequencer (ONT) operated using MinION Release v18.05.5 for the S1 run and 19.05.0 for the S2 run. The MinION sequencing device was connected to a computer running the linux operating system Ubuntu 16.04 LTS.

Data analysis

We used guppy for basecalling with the flip-flop model, trimming and homopolymer correction (guppy-cpu-v2.3.1). We then used guppy_barcoder for demultiplexing. For variant calling, we selected reads longer than 2000 bp with a mean quality higher than Q10 using NanoFilt [61]. The filtered sequences were then aligned using minimap2 on the Equus caballus genome reference assembly (GCF_002863925.1_EquCab3.0) [62]. We substituted the mitochondrial genome present in this assembly by a duplicated, pseudo-circularized, mitochondrial genome (GenBank JN398434) to increase the recovery of mapping to the circular mitochondrial genome. We used bioawk and samtools to extract mitochondrial reads and compute their quality scores [63]. These mtDNA reads were aligned with minimap2 on an Arabian breed reference horse mtDNA (JN398434), a Thoroughbred reference mtDNA (KC202971) or on the mtDNA from the reference Thoroughbred horse genome Twilight (MH586816). We selected medaka over nanopolish as the pipeline to call variant because our first tests showed that it required less computation time, and finally filtered the SNPs using a quality cut-off of 15. All visualisation scripts were written using the ggplot2 package in R3.6.0 [64, 65].

Sanger sequencing for SNP validation

We designed nine primer pairs to amplify nine fragments evenly spaced along the horse mtDNA (see Table 1). PCRs were performed on 4 ng of the same DNA used for long-read sequencing and using the 2X HotStartTaq Master Mix (Qiagen). PCR products were checked by agarose gel electrophoresis and purified using Nucleospin Gel and PCR clean-up kit (Macherey-Nagel). Cycle sequencing reactions were performed using the BigDye 3.1 chemistry, purified by precipitation and sequenced using an ABI 3130 sequencer (Applied Biosystems). Base calling was performed using the Sequencing analysis software version 5. 3 (Applied Biosystems).

The assessment of true positives (tp), false positives (fp) and false negatives (fn) was performed by inspecting the results of Sanger sequencing (on both strands) and of variant calling using the Geneious 11.1.5 software (https://www.geneious.com). We used a strict rule for our concordance analysis: we declared a SNP as a false negative or false positive if a difference was observed in one or more of the samples sequenced for a given horse, and even if there was a concordance in one or another sample. For each sequence run and each horse’s DNA we counted the occurences of tp, fp and fn. We then computed the values for the precision p = tp/(tp + fp) and the recall r = tp/(tp + fn). The F1 score was computed as the harmonic mean of precision and recall:

$$ \mathrm{F}1=\kern0.5em 2\ast \left[\frac{precision\ast recall}{precision+ recall}\right] $$