Oatk: a de novo assembly tool for complex plant organelle genomes

1Plastids and mitochondria are integral components of plant cells, governing photosynthesis and respiration, two fundamental physiological processes vital for plant growth and development. They both possess genomes of endosymbiotic origin that have undergone significant reduction and restructuring over evolutionary time. This includes gene losses, horizontal gene transfers, intracellular gene transfers, and inter- and intramolecular recombination, reflecting a complex evolutionary process shaping their current genomic structure and function. Plastid genomes (plastomes) of land plants are typically 120–160 kb in size and form a highly conserved circular, quadripartite structure, with a large single-copy (LSC) region and a small single-copy (SSC) region separated by two identical (or nearly identical) inverted repeats (IRs). In contrast, plant mitochondrial genomes (mitogenomes) demonstrate more extensive size variability, spanning a 200-fold range from tens of kilobases to over 10 megabases. They can exist in either circular or linear forms, often also contain long exact repeats, and may consist of multiple molecular components, highlighting the diversity and complexity of plant mitogenome architecture. Beyond their direct importance in plant physiology, organelle genomes provide valuable insights into genetic diversity and evolutionary relationships among plant species. Organelle DNA markers are widely employed in phylogenetic and evolutionary analyses and DNA barcoding initiatives to enhance our understanding of plant biodiversity and taxonomy.

2Organelle genomes have been sequenced using various approaches, starting from Sanger sequencing using primer walking, to sequencing from purified organelle DNA, and sequencing from whole-cell DNA. While the first two approaches are labor-intensive, costly, and sometimes species-specific, whole genome sequencing (WGS) provides a more efficient, robust, and cost-effective way to reconstruct plant organelle genomes. However, the long exact repeats in almost all plant organelle genomes, along with other issues discussed below, create severe difficulties for standard genome assembly tools. For this reason, there has been a requirement for dedicated software. Most of the contemporary tools for organelle genome assemblies are tailored for WGS data, such as MITObim, IOGA, NOVOPlasty, Organelle_PBA, and more recently GetOrganelle, GSAT, MitoHiFi, and PMAT. While the specific implementations may vary, two fundamental components persist in organelle genome assembly from WGS data: (1) distinguishing between organelle-derived and nuclear-derived sequences (potentially from old organelle genome integrations) and (2) assembling reads into genome sequences. The sequence differentiation usually involves the use of a seed database, which may consist of complete organelle reference genomes or sequence fragments derived from conserved genomic regions. Raw reads or assembled sequences are classified as organelle sequences if they overlap with these seed sequences. Given the potential divergence between seed sequences and the target genome, an “extension” step is often necessary. This step extends the sequences from the conserved regions with seed hits into neighboring divergent regions without seed hits by using sequence overlaps. The read assembly often starts with tools originally designed for nuclear genome assemblies, such as SOAPdenovo2 and SPAdes used by MITObim, IOGA, GetOrganelle and GSAT for assembling short reads, and the Newbler assembler employed by PMAT for assembling long reads.

3Despite the plethora of tools available for organelle genome assemblies, each has inherent limitations. First, most current tools rely on computationally intensive genome assemblers. Although strategies such as assembling only seeded organelle reads or only a limited subset of the WGS data are employed, the computational challenges remain substantial. Second, the “seed-extend” strategy can be problematic, particularly in assembling mitogenomes. Mitogenomes exhibit greater diversity compared to plastomes, posing challenges in creating a universal seed database that accurately represents a wide range of species. For example, GetOrganelle includes only four species for plant mitogenome seed database, leading to very limited representation. Notably, most existing tools, except for GSAT and PMAT, are primarily designed for plastomes and have limited application in mitogenome assembly, tested only on a restricted number of species. Third, many existing tools are optimized for short reads and lack compatibility with long reads, thereby missing a powerful resource for spanning repeats to elucidate complex genome structures. This limitation becomes more pronounced in light of the widespread adoption of highly accurate Pacific Biosciences circular consensus sequence (CCS or HiFi) in genome assemblies. Several tools have explored the use of long reads, but none fully capitalizes on the potential of PacBio HiFi reads: Organelle_PBA targets long noisy PacBio CLR reads, GSAT utilizes long reads only for resolving assembly graphs constructed from short reads, MitoHiFi is only aimed at simple circular mitochondria, and PMAT employs the Newbler assembler, which is slow and less suitable for high coverage data. The limitations of existing tools highlight the necessity for a dedicated organelle genome assembly tool that can effectively address these challenges. This tool should provide a more efficient and rapid approach to assemble organelle genomes from long-read sequencing data, as well as a robust method for accurate identification of organelle genome sequences from the whole-genome data to mitigate the misclassification and incomplete recovery problems, especially in the context of complex mitogenomes.

4Here, we present Oatk to address these challenges. Oatk is a de novo organelle genome assembly toolkit for assembling plastid and mitochondrial genomes from WGS data of high-accuracy long reads, specifically the reads used by many recent genome sequencing projects. It features a modular design, high speed, and user-friendly interface. Specifically, (1) we have developed a highly efficient genome assembler based on a sparse k-mer graph to mitigate the intensive computational demands for read assembly; (2) we have constructed hidden Markov model (HMM) profile gene databases to cover the full breadth of land plant organelle genomes and employ them instead of seed sequences for organelle sequence identification, achieving more accurate classification, particularly for mitogenomes; and (3) we have implemented a sophisticated graph resolution algorithm to generate primary assemblies, taking into account the graph structure and repeat-spanning sequence copy numbers. While the three modules can operate collectively through a wrapper program for organelle genome assemblies, each module can also function independently as a standalone command-line tool, offering enhanced flexibility: the genome assembler can produce assemblies for nuclear genomes comparable to MBG; the gene database facilitates gene annotation akin to MitoZ; and the graph resolution module can be utilized in a manner similar to GetOrganelle. We used Oatk to generate organelle genome assemblies for 195 species and compared its performance in read assembly with MBG and PMAT, and in graph resolution with GetOrganelle. Our findings suggest that Oatk generally outperforms the other tools across various metrics. Furthermore, we conducted a survey of the characteristics of the assembled organelle genomes and observed substantial genome diversity both within individual species and between species, for both plastid and mitochondrial genomes.

5Oatk consists of three components: (1) syncasm for genome assembly using a sparse de Bruijn graph, (2) hmmannot for sequence annotation based on a profile HMM gene database, and (3) pathfinder for graph resolution leveraging the graph structure and sequence coverage. Below we provide an overview of the key steps and concepts, with details given in the Methods section.

6For genome assembly, we first collect closed syncmers from the HiFi sequence data (Fig. 1a) and use them to build a sparse de Bruijn graph (Fig. 1b). Syncmers, invented by Edgar, provide a method for selecting a spanning set of overlapping representative k-mers, based on the position within them of their lowest scoring s-length substring (where ). Here we set default and. Next, the HiFi reads are mapped to the graph so as to identify low-frequency syncmers derived from sequencing errors and correct them. The graph is then reconstructed from the error-corrected syncmers, followed by graph cleaning and disentangling, and generation of a final assembly graph of unitigs (Fig. 1c). This unitig graph is a general assembly graph including the full nuclear genome. For organelle assembly, we next apply a k-mer coverage filter removing unitigs below five times the estimated nuclear haploid coverage depth to remove almost all nuclear sequence, including all but the most recent nuclear mitochondrial or nuclear plastid integration sequences (NUMTs/NUPTs) because k is large.

7For sequence annotation, each unitig sequence from the assembly graph is searched against a new land plants organelle gene database, comprising 130 plastid and 81 mitochondrial HMM gene profiles (Additional file 1: Tables S1 and S2). The resulting table of gene hits (Fig. 1d) is used to classify each remaining graph component into designated organelle types or non-organelle units (Fig. 1e).

8At this point, some of the graph components may still not be simple circles or loops, as standardly reported as organelle genome components, because of repeated sequences that can often be traversed in a variety of alternative orders (Fig. 1e). To reduce them to canonical form, we find the circular or linear path (or set of paths) that is best supported by spanning reads (Fig. 1f). The process to resolve this uses both the graph structure and sequence coverage, involving two major steps: estimating the copy number of each join (edge in the graph) between unitigs (vertices in the graph), and an exhaustive search through graph paths so as to maximize consistency with the edge copy number while covering all sequence in the graph (details in Methods section). This provides a minimal set of best supported structures that represents all the relevant sequence in the reads. We discuss below that in many cases there is evidence for alternate structures generated by recombination between repeats.

9We generated organelle genome assemblies for 195 land plants sequenced by the Tree of Life program at the Sanger Institute, mostly from the Darwin Tree of Life (DToL) project, including 24 monocots, 154 eudicots, 16 mosses, and one liverwort (Additional file 1: Table S3). In addition to Oatk, we ran MBG and PMAT for genome assembly construction to compare with the Oatk genome assembler syncasm. We also ran GetOrganelle for graph resolution to compare with the Oatk graph resolver pathfinder (see Methods section). The results are summarized in Table 1. Additional details can be found in Additional file 1: Table S4.

10Oatk generated both plastome and mitogenome sequences for all 195 species. The assembly graphs for these species are presented in Additional file 2: Fig. S1 for plastomes and in Additional file 2: Fig. S2 for mitogenomes. Oatk generated single-circular plastomes for all 195 species. For mitogenomes, Oatk generated 136 single-circular, 31 multi-circular, and 28 non-circular assemblies. Here we classified an assembly as multi-circular if it consists of multiple circular components, and as non-circular if it contains at least one linear component, regardless of whether it has a single component or multiple components.

11The results for assembly graph construction using syncasm and MBG are very similar, with the MBG assemblies resulting in just two fewer circular plastomes and six fewer circular mitogenomes than syncasm when using pathfinder for graph resolution (Table 1). At the base pair level, out of 195 species, 172 plastome and 131 mitogenome assemblies are identical, while 187 plastome and 157 mitogenome assemblies differ by no more than 3 bp (Additional file 1: Table S4). GetOrganelle consistently performed better with syncasm graphs than with MBG graphs, yielding more circularized assemblies, fewer non-circularized assemblies, and fewer failed cases. This suggests that syncasm graphs are more correct than MBG graphs in certain cases. This is anticipated, as syncasm has been optimized for organelle genome assembly, particularly for disentangling mixed organelle genome assembly graphs resulting from shared sequences (see Methods section, Additional file 2: Fig. S3). While syncasm and MBG delivered genome assemblies for all species, PMAT failed for 35 species due to either runtime errors or reaching the 24-h wall time limit. PMAT also shows a high error rate in organelle sequence classification. Out of the 160 species successfully assembled, we observed misclassifications in 144 species (Additional file 1: Table S5). This is likely because PMAT was initially designed for assembling mitogenomes. We therefore disregarded the classification results of PMAT and used only the assembly graphs as inputs to run pathfinder and GetOrganelle. Regarding speed performance, syncasm is slightly faster than MBG on average but has higher memory requirements (Additional file 1: Table S4). Both syncasm and MBG are significantly faster than PMAT, despite the extensive data downsampling for running PMAT (see Methods section, Additional file 1: Table S4).

12For assembly graph resolution, pathfinder consistently outperforms GetOrganelle irrespective of organelle types and genome assemblers, as reflected in the summary statistics. While the performance difference is marginal for plastomes, it is notable for more complex mitogenomes (Table 1). Misidentification of mitochondrial sequences is frequently observed with GetOrganelle, involving both the misidentification of plastid sequences as mitochondrial and mitochondrial sequences as non-organellar (Additional file 1: Table S4). For instance, GetOrganelle consistently misidentified the IR sequence on the plastomes of moss species as mitochondrial. This is probably because GetOrganelle relies on a label database with limited representation for plant mitogenomes, which may not adequately cover certain species. GetOrganelle was generally faster than pathfinder; however, the increased time is not a barrier to use of pathfinder. In most species, GetOrganelle completed within a few seconds, whereas pathfinder took up to a minute (Additional file 1: Table S4). It should be noted that, for pathfinder, the time for running hmmannot was included to ensure a fair comparison with GetOrganelle. Pathfinder alone typically completed within a second and successfully handled all assembly graphs that GetOrganelle failed to process (Table 1).

13The majority (182/195) of the assembled plastomes represent a standard quadripartite LSC-IRa-SSC-IRb structure which is clearly revealed by the assembly graphs (Fig. 2a–c). As described in Methods section, Oatk outputs these with the standard relative ordering of the large and small single copy regions (LSC and SSC respectively) as in the original Nicotiniana tabacum genome. The other 13 species demonstrate a simple circular assembly graph structure due to the absence of IRs (Fig. 2d), which all fall within previously identified inverted repeat lacking groups, including Erodium maritimum from the Erodium genus and twelve species from the legume identified inverted repeat lacking clade (IRLC). However, since these species are from several genera, their genome sizes vary considerably, ranging from 118 kb in Erodium maritimum to 149 kb in Lathyrus aphaca. We also observed a broad range of genome sizes for the species with quadripartite structures, from 122 kb in Lunularia cruciata to 233 kb in Schoenoplectus lacustris (Fig. 2e). Among them, bryophytes have relatively small genomes: the 13 mosses from the class Bryopsida have sizes ranging from 123 to 126 kb; the three mosses from the class Sphagnopsida have slightly larger genomes of about 140 kb; the sole liverwort, Lunularia cruciata, possesses one of the smallest genomes among all assembled species, at approximately 122 kb, second only to Erodium maritimum. In comparison to bryophytes, grasses (Poales) generally have significantly larger genomes, with 13 out of the 19 assembled genomes larger than 180 kb, including 12 sedges (Cyperaceae) and one rush (Juncaceae). The six species smaller than 180 kb are three rushes (Juncus squarrosus: 163 kb, Juncus bufonius: 170 kb, and Juncus effusus: 175 kb) and three grasses (Poaceae; Holcus mollis: 135 kb, Bromus sterilis: 137 kb, and Phragmites australis: 138 kb). Most of the other species have genome sizes lying between those of bryophytes and Poales, with 136 out of 146 falling within the range of 145 to 170 kb, four below this range, and six above it. Specifically, Calluna vulgaris and Jasione montana have genomes exceeding 200 kb, and mistletoe Viscum album represents the smallest genome among all the angiosperms assembled in this study, approximately 129 kb. The expansion of genome size generally resulted in an increased number of protein-coding genes, while the number of unique genes remained relatively stable, with a mean of 80.6 and a standard deviation of 2.0. Bryophyte species, particularly three Sphagnopsida mosses, are exceptions, with 90, 91, and 91 unique genes respectively. In contrast, Viscum album lost almost all genes related to NADH dehydrogenase, resulting in a unique gene count of 70, which is well below the average (Fig. 2e, Additional file 1: Table S6).

14Despite having the same quadripartite structure, significant variations in the sizes of different components were observed across species (Fig. 2a–c). Generally, genome sizes positively correlate with the sizes of the LSCs and IRs, but not with the SSCs. This is particularly evident in Poales, where we see notably expanded LSCs and IRs, accompanied by contracted SSCs (Fig. 2e). The average SSC size of the 13 Poales species with genomes exceeding 180 kb is 9.8 kb. In contrast, the average SSC size of the 157 species with genomes smaller than 180 kb, excluding Poales, is 17.9 kb. While this may appear counterintuitive, further investigations suggest that in larger plastomes the SSCs are partially duplicated and transferred into IRs. This IR boundary shift phenomenon is not exclusive to Poales species; it was also observed in other species including Jasione montana, Calluna vulgaris, Vaccinium vitis-idaea, Inga laurina, Inga leiocalycina, and Inga oerstediana, all of which have plastomes exceeding 175 kb. Moreover, the shifting of IR boundaries in these species is asymmetrical with respect to the two IR/SSC boundaries. The boundary adjacent to the ndhF gene remains relatively stable, with duplication always occurring from the opposite end of the SSC (Additional file 2: Fig. S4). The proportion of the SSC that was duplicated varies across species. In extreme cases, the SSC contracted to such an extent that only the ndhF gene remains, as observed in Jasione montana (4.1 kb), Vaccinium vitis-idaea (3.1 kb), and Calluna vulgaris (2.8 kb, Fig. 2c and g). The dynamics of plastomes extend beyond IR/SSC boundary shifts, encompassing a broad spectrum of phenomena. For example, in Berberis vulgaris, the IRs shift into the LSC, while in Luzula sylvatica, the IRs shift into both the LSC and SSC. In Hypericum pulchrum and Geranium rotundifolium, the LSC and SSC shift into the IRs, accompanied by intricate genome rearrangements, resulting in anomalous component sizes. In Viscum album, the SSC is highly contracted, leading to the loss of all genes related to NADH dehydrogenase (Fig. 2e, Additional file 2: Figs. S5–S9).

15In addition to the structural variation observed across different species, we also detected abundant structural heteroplasmy within individual species. From the point of view of the assembly graph, the single copy regions SSC and LSC can in principle be traversed in two distinct relative orientations, indicating that the plastome graph represents two distinct structural configurations (Fig. 1f). It has long been reported that both forms can coexist within a single individual, a phenomenon that can be explained by flip-flop recombination between two IRs. Recently, Wang and Lanfear revisited this discovery using a method combining long-read mapping and statistical modeling and confirmed that this type of heteroplasmy exists in most plants. Here we used assembly graphs to facilitate a more direct investigation of this topic. Since a read must map through the entire IR region in the graph to provide evidence for a specific form, we set the IR size limit to 20 kb taking into account the length restriction of PacBio HiFi reads, resulting in a selection of 22 species, including 17 bryophytes and five angiosperms. For each species, we mapped the HiFi reads to the individual assembly graph and count the number of reads supporting each form. For all species tested, we see approximately the same number of reads supporting both forms (Fig. 2f) with no significant deviation from a binomial distribution with probability 0.5 (Additional file 1: Table S7). This observation reaffirms previous findings that the two forms widely coexist in plants, with their ratio being close to 50/50 and is consistent with frequent intramolecular flip-flop recombination as a plausible explanation for this phenomenon. Besides chromosomal-scale structural variations, we also observed other types of heteroplasmy at smaller scales. Figure 2c shows the assembly graph of Calluna vulgaris as an example. The IR region contains two nested small repeats, leading to a complex structure. By mapping the reads to this region, we identified multiple isoforms with varying molecule abundances (Fig. 2h). Similar heteroplasmic structures were observed in several other assembled species (Additional file 2: Figs. S1, S10). However, it is challenging to distinguish between non-identical IRs within a single molecule and heteroplasmy between molecules. To confirm non-identical IRs within a molecule, a read must span a single-copy region and cover both variants, which is difficult to observe with HiFi reads, as they are typically shorter than 20 kb, but should be possible with longer ONT reads.

16The assembled mitogenomes are much more variable than the plastomes, demonstrating remarkable diversity in size, structure, and gene content (Fig. 3a–g, Additional file 2: Fig. S2). Genome sizes range from 103 to 2731 kb, with the smallest genome observed in the moss Thuidium tamariscinum and the largest in the sedge Carex laevigata. A similar trend in genome sizes was observed in mitogenomes as in plastomes, with bryophytes having relatively small mitogenomes, Poales mitogenomes being notably larger on average, and other species falling between them. There is a positive correlation between gene number and mitogenome size in general (Fig. 3g, Additional file 1: Table S8). Among the 195 mitogenomes analyzed, 136 represent a single circular structure. These single circular mitogenomes exhibit various underlying structures, including simple ring structures, quadripartite structures with either an inverted or directed repeat (DR), and more complex structures involving multiple circular molecules linked by repeat sequences (Fig. 3a–d). Thirty-one of the assembled mitogenomes consist of multiple circular components (Fig. 3e). Notably, the Galeopsis tetrahit mitogenome exemplifies this complexity, comprising 368.8 kb of sequence distributed across ten circular molecules, each containing protein-coding genes (Additional file 2: Fig. S11). The remaining 28 assembled mitogenomes are non-circular, comprising either a single linear component or multiple components, with at least one being linear. While some of these non-circular mitogenomes are probably the consequence of incomplete assembly, others can be verified through read mapping. For instance, the Alnus glutinosa mitogenome consists of two components: one circular and one linear, with the latter displaying an intriguing loop structure at each end (Fig. 3f). Reads were found mapping through each direction of the bifurcation in the loop structure, with a smooth change in coverage, confirming the authenticity of this structure (Fig. 3h). The loop structure is not unique to Alnus glutinosa; similar structures were observed in at least seven other species assembled in this study (Additional file 2: Fig. S2).

17Structural heteroplasmy was commonly observed in the assembled mitogenomes. The genome of Pseudognaphalium luteoalbum exhibits a standard plastid-like quadripartite structure with IRs (Fig. 3b). Similar to observations in the plastome, this structure manifests in two coexisting forms, characterized by different orientations of the two single-copy sequences. The number of reads mapping across the IRs supporting the two forms are 299 and 301, respectively, suggesting unbiased interchange. We observed a similar phenomenon in the genome of Calluna vulgaris, which also represents a quadripartite structure but with DRs (Fig. 3b). The number of reads supporting the two forms are 328 and 361, respectively. Unlike the IR-featured structure, in this scenario, the second form consists of two separate small circles, supported by 170 and 191 reads, respectively. We extended this analysis to all accessible bidirectionally bifurcated repeats across the assembled mitogenomes. Specifically, we examined each sequence on the assembly graphs and included it as a repeat in the analysis set if it meets the following criteria: (1) the sequence is no larger than 10 kb; (2) the sequence has two incoming edges and two outgoing edges; and (3) the size of each sequence associated with the four edges is at least 10 kb. The sequence size threshold was set to ensure the quality of read mapping. This process identified 166 repeats across 109 genomes. It is important to note that, because we use a k-mer size of 1001, all these repeats are larger than 1 kb. The read mapping results suggest that nearly all these repeats are associated with structural heteroplasmy, with many representing nearly balanced abundances (Fig. 3i, Additional file 1: Table S9). The repeat-driven dynamics of local structures can extend globally, leading to diverse molecular forms involving various combinations of subgenomic regions (Fig. 3j). The coexistence of these superstructures and subgenomic structures has been well documented for decades and characterized using various analytical techniques, with one possible explanation being intramolecular recombination.

18We observed frequent sequence sharing between plastid and mitochondrial genomes within a species, reflecting historical transfer of DNA from one organelle to the other. This usually led to a mixture of sequences from two organelle genomes in the same assembly graph forming a entangled structure. Figure 4a demonstrates the genome assembly graph of the Alnus glutinosa as an example. It contains four identical shared sequences with sizes from 1382 to 4753 bp. The homologous sequence pair forms a bubble structure, differentiated by a few variants, including seven SNPs and three 4–7 bp small INDELs. For each species, we aligned the plastid and mitochondrial genomes to identify shared sequences. We set a minimum sequence size of 1 kb and a minimum sequence identity of 90% to eliminate noisy alignments. We identified 578 shared sequences totaling 1.57 Mb across 144 species, with individual shared sequence sizes ranging from 1 to 17.3 kb. Among them, 59 species contain more than 10 kb of shared sequences, with Carex laevigata having the most at 66.6 kb. This finding demonstrates the prevalence of sequence transfer between organelle genomes, although, notably, no shared sequences are found in any of the bryophyte species. No clear correlation between mitogenome size and the amount of shared sequences was observed (Fig. 4b). For example, Erodium maritimum has a 307.7 kb mitogenome and 37.9 kb of shared sequences, whereas Juncus bufonius has a 993.7 kb mitogenome but only 3.4 kb of shared sequences. The majority of sequence transfers occurred from plastome into mitogenome, not the other way around. Out of the 1.57 Mb of shared sequences identified, 1.11 Mb were annotated as gene sequences, with 1.09 Mb (98.4%) specifically identified as plastid genes (Additional file 1: Table S10).

19In terms of sequence similarity between paralogs, 288 sequences totalling 902.3 kb show greater than 98% identity, while 106 sequences totaling 222.2 kb show less than 95% identity. Although many paralogs represent high similarities, sequence identities are broadly distributed across the entire range from 90% to 100%, indicating that sequence transfer is a repeated, ongoing process (Fig. 4c). This is supported by phylogenetic analysis of 65 paralogous sequence pairs annotated as the psaA gene—the most frequently annotated shared gene followed by psbB with 60 sequence pairs. The phylogenetic tree constructed from these 130 sequences accurately represents the clade structure and shows a stratification of shared sequences (Fig. 4d). For many paralog pairs, two sequences were closely clustered together in the tree. Notably, sequences in 19 paralog pairs are mutually the nearest neighbor of each other and sequences in 33 paralog pairs are mutually within the three-nearest neighbors (Additional file 1: Table S11). Clades for some genera like Hypericum and Inga show separate clusters for plastid and mitochondrial sequences, suggesting that sequence transfer occurred before their speciation. In 50 out of the 65 paralog pairs, the mitochondrial sequences demonstrate greater distance to the node representing the most recent common ancestor (MRCA) compared to their plastid paralogs, such as observed in species Juncus squarrosus, Chamaenerion angustifolium, and Iris foetidissima (Additional file 1: Table S12). Statistical analysis using the Wilcoxon signed-rank test on all 65 paralog pairs indicates a significant difference in distances (p value = ). This presumably reflects loss of selective constraint in the mitogenome.

20We developed Oatk for de novo assembly of plant organelle genomes using high-accuracy long reads and applied it to assemble 195 species spanning a broad range across the tree of life. Compared to other tools, Oatk produced complete assemblies for more species, particularly for mitogenomes. Oatk is distinguished from other tools by a highly efficient genome assembler utilizing a sparse k-mer graph to reduce memory usage and accelerate graph construction. For k-mer selection, we use closed syncmers, which are intrinsically suitable for k-mer sparsification. The sparsification factor is approximately, where k and s are the k-mer and s-mer size respectively. With high-accuracy long reads, we can choose a large k and typically have. In this case, the average distance between adjacent k-mers is roughly k/2. For example, with as used in this study, the data is compressed by a factor of about 500, thereby facilitating efficient assembly graph construction. Other techniques for k-mer sparsfication include the minimizer approach used in MBG and the k-min-mer approach used in MDBG. MBG features a syncmer-based implementation when integrated into Verkko. The major difference between syncasm and the syncmer-based MBG probably lies in graph disentangling, with syncasm being better tailored for organelle genome assemblies.

21Even though there are many well-designed genome assembly tools for HiFi data, such as hifiasm and HiCanu, they have not been effectively used for organelle genome assemblies. These assemblers are not optimized for extremely high coverage data, such as organelle genomes, where coverage can reach thousands or tens of thousands of times, leading to problems with sequence error correction. Additionally, the highly unbalanced data coverage between organelle genomes and nuclear genomes presents another challenge. Finally, these tools generally rely on string graphs, which suffer from inaccurate sequence coverage estimation—a critical factor for graph resolution in plant organelle genome assemblies. In contrast, k-mer based assembly graphs provide better sequence coverage estimates. Given a long read assembly with coverage data, GetOrganelle takes a similar approach to Oatk for finding the optimal path, but Oatk deals with some additional cases that GetOrganelle fails on; this appears to be in part because Oatk has a more extensive labeling database so correctly identifies more graph components, and in part because in its second phase path search Oatk adjusts edge coverages rather than node coverages, leveraging some additional information.

22Another distinctive characteristic of Oatk is its use of a profile HMM gene database for organelle sequence identification, rather than relying on seed sequences. For a given gene, a profile HMM was constructed considering all related sequences available in the NCBI repository to ensure a broad species representation. This approach provides a more sensitive method for sequence identification, particularly useful for divergent mitogenomes. Oatk is not the first tool to use profile HMMs for genome annotations. They have been employed by MitoZ to identify target sequences in animal mitogenome assembly, and by GeSeq as an essential component for protein and rRNA coding gene annotation in plant organelle genomes. Although our focus has been on plant organelle genome assemblies, Oatk extends its utility to other species with the corresponding gene database. We have created HMM profile gene databases for various clades and developed a tool to facilitate the creation of gene databases given an NCBI taxonomy ID (Methods section). As an example, we applied Oatk to assemble complete mitogenomes for five animal species, including a mammal, a fish, a bird, a lizard, and an insect (Additional file 2: Fig. S12). This demonstrates the broader applicability of Oatk beyond plants. OatkDB maintains HMM databases for a wide range of lineages and offers user-friendly scripts for constructing custom databases (Methods section). Nevertheless, it is important to acknowledge that not all organelle genes will necessarily be included, particularly for lineages with limited well-assembled organelle genomes. A potential future direction could be to incorporate a de novo gene annotator.

23While Oatk has demonstrated high efficiency and reliability, it has some limitations that should be considered. Oatk uses a very large k-mer size, making it optimized for high-accuracy long-read data such as HiFi reads, ONT duplex reads, or error-corrected ONT reads. Since its k-mer size typically exceeds the length of short reads (e.g., Illumina), Oatk is not suitable for short-read data. Although theoretically possible with a smaller k-mer size, this would likely result in an intractable assembly graph, entangling organelle and nuclear genomes due to the nonuniqueness of small k-mers. GetOrganelle mitigates this issue by identifying organelle sequences at the sequence level and limiting the assembly process to predominantly organelle sequences. Oatk is also less suitable for noisy long-read data, such as those generated by PacBio CLR or uncorrected ONT technologies, as their typical error rate would lead to an unacceptable level of erroneous k-mers when k is large. However, error-correction tools such as HERRO and hifiasm can process ONT R10 reads, allowing error-corrected ONT reads to be used as input for Oatk. We note that another recently published method TIPPo is also limited to long high-accuracy reads. Another limitation is that Oatk employs coverage-based pre-filtering to remove nuclear genome sequences, which relies on adequate excess coverage in organelle genomes. We allow the user to alter the default threshold of five times average nuclear coverage, but Oatk would not be able to handle cases where organelle coverage is comparable to or less than nuclear coverage. Finally, Oatk does not explicitly assess the completeness of an assembled organelle genome. However, it provides the assembly graph and gene annotation results, which can assist users in making this determination. For example, a circular genome assembly containing all core genes is more likely to be complete.

24We observed extensive diversity in the assembled genomes, both within and across species. Analyses combining long reads and assembly graphs offered significant insights into the structural heteroplasmy of these genomes, highlighting their potential as powerful tools for understanding complex organelle genome structure. However, this poses an important question regarding the appropriate representation of plant organelle genomes. Traditionally, genomes are stored as linear sequences and labeled as circular if applicable, as required by public sequence repositories in the International Nucleotide Sequence Database Consortium (INSDC: NCBI GenBank, ENA, and DDBJ), whereas the structural heteroplasmy information present in organelle genomes is overlooked. Some tools, such as GetOrganelle and GSAT, opt to output sequences of all possible conformations, a practice Oatk could potentially adopt. However, this approach introduces additional challenges such as sequence redundancy and combinatorial complexity, especially in the case of complex mitogenomes where the number of conformations may escalate to an unmanageable level. We suggest that it might be more useful to provide a primary assembly that includes all sequences together with the assembly graph, allowing flexibility depending on the research question. This underscores the need for tools that can directly analyze assembly graphs, such as for gene annotation, where considering graph structure may change annotation outcomes.

25Oatk provides a new, robust approach to the assembly of plant organelle genomes, suitable for use by large scale reference genome sequencing projects. For example, it is now standardly used by the Darwin Tree of Life Project for plant species assemblies. In addition to outputting an optimized supported circular genome, or where necessary a set of circular and linear genome components, Oatk also provides a full genome sequence graph that represents all supported heteroplasmic variants.

26Closed syncmers are parameterized by two positive integers s and k and a hash function of strings of length s that returns the same value for a string and its reverse complement. Given a string of length k, also called a k-mer, to determine if is a closed syncmer, is first used to compute the hashes for all s-long substrings (s-mers) from. The k-mer is a closed syncmer if and only if the hash of the first or last s-mer of is minimal, i.e., it is not larger than all other hashes from. We use and in this study.

27Given a set R of HiFi sequencing reads, we build a k-mer table T containing all the distinct closed syncmers present at least once in the reads. Since most HiFi sequencing errors are homopolymer run length errors, the reads are homopolymer compressed (hoco) before collecting syncmers, where homopolymer runs are collapsed into a single nucleotide. In T, each k-mer is assigned a unique sequential index. The table also maintains the the k-mer frequency (coverage), its hoco sequence, and the its full length consensus sequence. To generate the consensus sequence of a hoco k-mer, the homopolymer run lengths for each base of the k-mer are recorded during the k-mer table construction. The run length consensus of a base is computed as the closest integer to the average of all run lengths at its respective positions, mirroring the approach used in MBG. We use a hash function, denoted by, to map k-mer sequences to 64-bit integers for faster comparisons during the table construction. It is similar to for s-mers, excepting that while is designed to be a perfect hash function eliminating hash collisions, hash collisions for are inevitable, owing to the values of s and k. The hash collisions are checked, and distinct table entries are created for k-mers with collided hashes. To allow for both directions of double-stranded DNA, is applied to compute hash values for both the original k-mer and its reverse complement, and only the sequence with the smaller value constitutes a k-mer entry in the table.

28Given the k-mer table T, each read in the read set R can be rewritten as a vector of syncmers accompanied by the respective positions on the read and orientations concerning the DNA strands. Let be a HiFi read harboring m closed syncmers, denoted by sequentially. The respective positions of these syncmers on r are in ascending order. Let be a function mapping a k-mer to its index in the table. Let be a function calculating the relative orientation of a k-mer compared to the corresponding entry in the table, with 0 being the same strand and 1 being the reverse complement. The read r can then be rewritten as, where is the transformation function. Let be the transformed read set after applying to each read in R. The inherent characteristic of assured overlap between consecutive closed syncmers, i.e., for all, guarantees that the whole read apart from material before the first syncmer and after the last syncmer can be restored from its respective syncmer vector.

29Constructing a sparse k-mer graph from closed syncmers

30The k-mer table T and the transformed read set Z are used to construct a sparse k-mer graph. For each k-mer entry in the table T, two vertices are established, one for the k-mer itself and the other for its reverse complement; denote these as v and respectively as an example. Vertices v and are called complement vertices to each other, and holds. Vertices maintain references to the k-mer indices, allowing for retrieval of k-mer information from the table, including details such as the sequence and the coverage. Two complement vertices refer to the same k-mer entry in the table. However, it is necessary to calculate the reverse complements when retrieving the underlying k-mer sequence and k-mer consensus sequence for the one not present in the table.

31Directed edges are introduced between vertex pairs when the corresponding k-mers are adjacent on a transformed read. For any edge, a complement edge is also added. Consider the example of the read. Let v and w be the associated vertices of the first two k-mers and, respectively. Given the read r, a directed edge will be added to the graph, pointing from vertex v to w, and simultaneously its complement edge is added, where is the respective complement vertex of v. Edges record their coverage representing the count of positions on the transformed reads supporting the adjacency of the two vertices. In the example, the read contributes coverage of one unit to both the edges and. In addition, edges retain the length of the sequence overlap between the k-mers represented by the two vertices. For a given edge, the sequence overlap is determined by consolidating all overlap distances calculated from the positions of the two corresponding adjacent k-mers on the transformed reads. Typically, all the overlaps for a particular k-mer pair edge will be identical, and this uniform distance will be utilized to calculate the overlap length. However, in rare cases, it is possible to observe multiple overlap lengths due to complex repeat structures; when this happens the largest distance will be employed to calculate the overlap length. The complement edges and exhibit identical edge coverage and sequence overlap length by definition.

32Let be the k-mer graph, where V is the vertex set and E is the edge set. The above construction procedure guarantees G is Watson-Crick complete: (i),; and (ii),, where and are the complement vertices of v and w, respectively. For all, let be the outdegree of v and be the indegree. It follows that.

33The k-mers in the table T with coverage below a predefined threshold are identified as potential error syncmers. Let be the set of potential error syncmers. Firstly, the vertices corresponding to the k-mers in are deleted in the k-mer graph G. Denote by G′ the trimmed k-mer graph, assumed to be free of syncmer errors. Then, for each read in Z, error blocks, defined as consecutive blocks of error syncmers, are identified and mapped to G′ for error correction. Consider the example of the read, let be an error block of b syncmers, where. By definition of the error block, (i); (ii) if; and (iii) if. For simplicity, we have omitted the mapping functions and here for converting actual k-mers to k-mer table entries. We call the left boundary of the error block and the right boundary. If both the left and right boundaries are absent, the entire read constitutes an error block and will be excluded from error correction. Without loss of generality, we assume the existence of the left boundary; otherwise, we can establish the left boundary by converting the read to its reverse complement. The syncmer vector delineates a subsequence of the read r harboring sequencing errors, starting from the position and concluding at the position. Let be the subsequence, and be the length in base pairs of the subsequence. In broad terms, the error correction aims to find an unambiguous syncmer path over G′ starting with the vertex corresponding to the left boundary outlining a sequence with a sufficiently small Levenshtein edit distance to. If the right boundary exists, the syncmer path must end with the vertex corresponding to it. The upper limit of the Levenshtein distance is parameterized by and calculated as for a subsequence of length l. The parameter should reflect the overall base accuracy of the sequencing data. In this study, we set and.

34To find the error-correcting syncmer path for an error block, we perform a depth-first search (DFS) on G′ starting from the vertex corresponding to the left boundary, denoted by. Let be the current search path of the DFS at depth n. Let be the corresponding sequence delineated by the search path, and be the sequence length. A dynamic programming (DP) table is first constructed for calculating the Levenshtein distance between and. Let be the edit distance of the prefix alignment between two sequences, i.e., the edit distance ignoring the trailing indels in the DP table. Let be the length of the overhang sequence of, i.e., the number of trailing insertions of in the DP table. The actual edit distance between and is computed as. The “depth search” along the path is concluded at under the following conditions: (i), or (ii), or (iii) there is no succeeding vertex from. When a “depth search” is completed, the subpaths along the search path are examined to update the global minimum and the second minimum edit distances, denoted by and, respectively. The path with the minimum edit distance is also updated. The DFS search routine then proceeds to the next search path. A constraint on the total number of “depth searches,” parameterized by, is set to prevent excessive path exploration in complex graph regions arising from highly repetitive sequences. The error correction is considered unsuccessful if the number of searched paths exceeds the constraint. Upon the successful completion of the DFS, the values of and are examined. The error correction is considered valid if and, where is a parameter to control path ambiguities, with lower values indicating higher confidence of the optimum path. For a valid error correction, the path with the minimum edit distance is used to update the syncmer vector of the read. In this study, we set and.

35The bottleneck in error correction arises from the computation of Levenshtein distance. We implemented two strategies to address this issue. Firstly, we employ Ukkonen’s O(ND) algorithm, where N is the minimum length of two sequences and D is the edit distance between them. The fundamental idea involves constraining the DP matrix computation to a band along the diagonal that is 2D wide. This implementation is intrinsically suitable for our needs, considering we set an upper limit of for the Levenshtein distance given a sequence of length l. Secondly, we reuse the DP matrix during the DFS. When visiting a vertex at each DFS step, the target sequence is extended by appending the sequence encoded in the current vertex, and the computation of the DP matrix is confined to only the part relevant to the newly added bases.

36Given that organelle genomes typically exhibit significantly greater coverage than the nuclear genome in the sequence data, we set a lower threshold for k-mer coverage, denoted by c, to exclude k-mers derived from the nuclear genome. After filtering the low-copy k-mers, we apply a strategy similar to Miniasm to remove tipping sequences shorter than ten kilobase pairs and pop bubbles shorter than 100 kilobase pairs. A unitig is defined by a “maximal unambiguous path” in the k-mer graph. More precisely, a vertex path forms a unitig if for and (i) or (ii) and. Since each vertex in the k-mer graph corresponds to a k-mer entry in the k-mer table, a unitig can be written as a syncmer vector, similar to a read. Let U be the set of unitigs. For, let be the syncmer vector for u, using the same notations for reads. A unitig graph is derived using U as vertices, and for, adding a directed edge pointing from to if the last k-mer of is identical to the first k-mer of. The structure and characteristics of the k-mer graph persist in the unitig graph. Specifically, the unitig graph remains Watson-Crick complete. To generate the sequence for a unitig, the k-mer sequences along the k-mer vector are concatenated, eliminating the overlapping sequences between each pair of adjacent k-mers. The final genome assembly comprises the sequences of all unitigs. For a complement unitig pair, only one copy is retained in the assembly to avoid redundancy.

37Horizontal transfers of DNA between organellar genomes can result in an entangled assembly graph for the two organelle genomes owing to the presence of shared sequences longer than the k-mer used for graph construction. To tackle this problem, we leverage the phase information of k-mers obtained from reads represented by vectors of k-mers. Firstly, the reads are mapped to the unitig graph at the k-mer level employing a seed-and-extension approach: the first k-mer of the read presented in the unitig graph is used as the seed to initiate a DFS; the DFS path is then extended based on the order of k-mers on the read. Then, for each unitig node, the number of reads supporting each spanning triplet are counted and used to resolve the local graph structure around the node. A spanning triplet of a node (the center node) is a graph path consisting of three nodes with an incoming and an outgoing node connecting the center node. For example, node u has incoming nodes, written as, and outgoing nodes, written as, then there are spanning triplets for u. Write, as the number of reads supporting the triplet path. The triplet path is called dominated if and, where is a predefined threshold and set as 0.1 in this study; otherwise, it is non-dominated. A unitig node is called resolvable if there exists at least one dominated spanning triplet. If a node is resolvable, we update the unitig graph by introducing a compound node for each non-dominated triplet. Precisely, for a non-dominated triplet, the compound node, denoted by, is constructed by concatenating the k-mer vectors of and, i.e.,; all the incoming edges pointing to will be copied and point to; all the outgoing edges pointing from will be copied and point from. The resolved node u is finally removed from the unitig graph. This resolution process is repeated for multiple rounds until the graph structure converges.

38We construct the profile HMM gene databases separately for plastid and mitochondrial genomes. Given the species and organelle genome type, we download all related sequences from the NCBI repository in Genbank format using the Entrez Direct command tools. Specifically, we use the esearch tool with the query string “TAXID [Organism] AND OTYPE [Filter]” to search against the nucleotide database for a list of target sequences and the efetch tool to download those sequences, with “TAXID” being “txid3193” for Embryophyta, “txid40674” for Mammalia, “txid7898” for Actinopterygii, “txid8782” for Aves, “txid8504” for Lepidosauria, and “txid50557” for Insecta, and with “OTYPE” being “chloroplast” or “mitochondrion.” A sequence set is created for each gene, including protein-coding genes, tRNA genes, and rRNA genes, after parsing the Genbank file and used to construct a profile HMM. For profile HMM construction, we first perform an initial sequence clean by removing sequences with invalid characters and those shorter than one-third or longer than three times the average length. If the number of sequences is more than 10,000 after the initial clean, a second-round clean is performed to retain the top 10,000 most representative sequences, determined by 12-mer completeness: a 12-mer table is first generated for all sequences, excluding the 1% most and 1% least frequent ones; a score is then computed for each sequence as the sum of counts for all unique 12-mers; the 10,000 sequences with the highest scores are finally chosen. With the cleaned sequences, the multiple sequence alignment is constructed using MAFFT (version 7.505) with the “--auto” option and the profile HMM is built using HMMER’s “hmmbuild” tool (version 3.4) with default parameters. All gene profile HMMs are concatenated into a single file and compressed with HMMER’s “hmmpress” tool (version 3.4) to make a binary database.

39The assembled sequences are searched against the profile HMM database with HMMER’s “nhmmscan” tool (version 3.4) to build a table for protein-coding genes. Gene hits with e values greater than or scores less than 300 are considered insignificant and omitted in the further analysis. The gene annotations are used to determine the organelle type for each assembly graph component. We construct two gene tables for each graph component, one for plastid genes and one for mitochondrial genes, containing only the best hit of each gene on the sequences belonging to the component. We could simply assign the organelle type of a component by comparing the sum of scores to plastid and mitochondrial genes. However, this sometimes leads to misassignments of mitochondrial sequences into plastid due to the frequent insertion of plastid sequences into mitochondrial DNA. Since plastid genes are well conserved, they have higher annotation scores than mitochondrial genes on average. As a result, a few inserted plastid genes can become dominant. To mitigate this problem, we use a progressive assignment approach rather than a simple comparison of the sum of scores. The graph components are sorted by sum of scores in descending order and processed progressively: (1) if the component has not been processed, assign it the corresponding organelle type; (2) if the component has been processed and classified as plastid, reassign it to mitochondrial if (i) there has already been a component assigned to plastid, and (ii) the mitochondrial score is no smaller than one-third of the plastid score. The above process could incorrectly identify some junk sequences as organelle sequences, including those originating from NUMT/NUPT or resulting from sequencing errors. To address this issue, we track the highest score for each gene in sequences already assigned to an organelle type. When we need to assign an organelle type to a new component, we compare the gene scores in sequences from this component to the current highest scores. A component is assigned only if it has a sufficient number of gene hits with higher scores than the current best. We set this threshold at three for plastid and one for mitochondrion. The highest scores for the genes are updated each time a new component is added.

40The genomes of plant organelles frequently cannot be assembled into a single contiguous sequence because of large repeats or alternative isoforms. We examine all possible path conformations within the assembly graph, selecting the longest path as the primary genome assembly. If multiple paths are of equal length, we choose the one most strongly supported by the read data, determined by the total number of reads spanning the node junctions along the path.

41We start by estimating the copy number of each sequence in the graph using the same method that GetOrganelle employs for calculating “observed multiplicity”. Consider an assembly graph of n nodes, denoted by,...,. The length of is and the sequence depth of coverage of is,. We estimate the average base sequence coverage, denoted by, and determine the copy number of as. This estimation is performed using an EM algorithm. In the E-Step, we calculate, and in the M-Step, we calculate. With the copy number information for each sequence, we aim to traverse the entire graph while ensuring each node is visited according to its copy number. To achieve this, we expand the graph by creating additional copies of multi-copy nodes and their corresponding edge connections. In the expanded graph, each node only needs to be visited once. We implemented the path exploration algorithm to ensure that: (1) all paths are examined; (2) equivalent paths due to rotation and reversal are examined only once; and (3) equivalent paths resulting from permutations of different copies of the same nodes are examined only once. During the path exploration, the circular and linear paths encompassing the longest sequences are added to the candidate set for the ultimate selection. For the final choice of the optimal path, a circular path is preferable over a linear path when the circular one covers no less than a predefined portion (90% in this study) of the sequences covered by the linear one in length; among the paths of the same length, the path with the most read support (designated as the graph edge coverage) is preferable over the others.

42In plastome assemblies, there are typically two possible path conformations due to IRs. We select the one that shows better gene order consistency with the Arabidopsis thaliana reference plastome, which has the same relative order as the original published Nicotiniana tabacum plastome. Specifically, we use the order of the 71 protein-coding genes from the A. thaliana assembly as a reference. We compute Spearman’s rank correlation coefficients between the reference and each of the path conformations and select the one with the higher correlation coefficient. Additionally, we rotate the assembled sequence so that it begins with the start codon of the psbA gene as a convention.

43If the primary sequence assembly resulting from the above process is non-circular, we adjust the sequence copy number estimates using the graph topology and undertake a second round of path exploration. Consider a node with incoming nodes, written as, and outgoing nodes, written as. Write and the number of visits to the respective edges connecting in the path, where and. Write and the real indegree and outdegree of the node. If the entire graph can be traversed with a Hamiltonian cycle, then for any node in the graph, holds. It is important to note the difference between and, while the former considers only the graph structure, the latter also considers the copy number of nodes. We minimize the objective function:. Intuitively, we minimize the sum of the differences between the indegrees and outdegrees, as well as the discrepancies between the expected node copy number and the observed degree. In the objective function, the term for each node is weighted by the logarithm of the corresponding sequence length to assign greater importance to longer sequences. The variables take discrete values and are bounded by. If the solution space is no greater than 100 million, we perform a brute force search; otherwise, we use simulated annealing.

44We ran Oatk to generate organelle genome assemblies for 195 land plant species sequenced by the Tree of Life program. For organellar sequence identification, the Embryophyta plastid and mitochondrial profile HMM gene databases were specified using the “-p” and “-m” parameters, respectively. GenomeScope (version 2.0) was employed to estimate the haplotype data coverage of the nuclear genome. The coverage threshold (“-c”) was set to five times the estimated haplotype coverage, and the number of threads (“-t”) was set to 6 for running Oatk, with all other parameters left at their default settings.

45We ran MBG (version 1.0.16) and PMAT (version 1.5.3) for genome assembly graph construction to compare with the Oatk assembler syncasm. For MBG, we used the same k-mer size and coverage threshold as syncasm, i.e., “-k” 1001 and “-c” five times the haplotype coverage, and set “-t” to 6 to use six threads. All other parameters were left at their default settings. For PMAT, we configured the genome size parameter (“-g”) to the estimated value from GenomeScope, set the sequence type parameter (“-st”) to “hifi,” and adjusted the sampling ratio factor parameter (“-fc”) based on the GenomeScope haplotype coverage estimation to use data. Additionally, we set the “-cpu” parameter to 6 to use six threads and enable the “-m” option to retain sequence data in memory for faster processing. PMAT generates separate assembly graphs for plastid and mitochondrial genomes. However, sequence misclassification often results in a significant proportion of shared sequences between the two assembly graphs. To address this issue, we merged the two graphs to create a unified assembly graph for both organelle genomes, removing redundant sequences and edges while preserving the overall graph structure.

46We ran GetOrganelle (version 1.7.7.1) to generate organelle genomes from assembly graphs constructed by syncasm, MBG, and PMAT. We set “-F” to “embplant_pt” and “embplant_mt” respectively for assembling plastid and mitochondrial genomes, and set “-t” to 6 to use six threads. We ran the Oatk hmmannot-pathfinder pipeline to compare with GetOrganelle. This was achieved by using the “-G” parameter in Oatk, which allows it to take an existing assembly graph as input. By specifying this parameter, Oatk bypasses the assembly graph construction from the raw sequence data.

47We ran NUCmer (version 4.0.0rc1) with the parameter “--maxmatch” to generate a self-alignment for the genome sequence. Alignments with sequence identity below 95% and those not in the reverse-complement orientation were filtered out. Alignments within a distance smaller than 50 bp were merged, where the distance between two alignments is defined as the Manhattan distance between their closest endpoints. Following merging, the two largest aligned segments are identified as IRs. The regions flanking these IRs are identified as LSC and SSC based on their respective sizes. The genome sequence is finally rotated so that it starts at the beginning of the psbA gene, which is at the start of the LSC, and ends with an IR.

48We ran GraphAligner (version 1.0.19) to map reads to the assembly graphs with the parameters “--precise-clipping 0.9 --min-alignment-score 5000 -x vg.” Alignments were filtered out if less than 95% of the read was aligned to a unique position. GraphAligner outputs a path on the graph for each alignment. For any given graph path of interest, we examine the alignment path of each read and count it as a supporting read if the alignment path fully covers the target graph path. Alignments in both directions are considered, and the sum is taken as the total number of supporting reads.

49We ran NUCmer (version 4.0.0rc1) with the parameter “--maxmatch” to map the plastome to the mitogenome. We then ran delta-filter (version 4.0.0rc1) with the “-l 1000” parameter to remove short alignments, the “-i 90” parameter to retain only alignments of at least 90% sequence identity, and the “-r” parameter to keep only the best mitogenome match to each region of plastome sequence. The sequences corresponding to the resulting alignments are considered the shared sequences between the two organelle genomes.

50We ran hmmannot to generate gene annotations of the shared sequences and extracted sequences annotated as psaA genes. We ran MAFFT (version 7.505) with the parameters “--localpair --maxiterate 1000” to generate a multiple alignment of these gene sequences and ran IQ-TREE (version 2.3.4) with default parameters to generate a phylogenetic tree. The calculations of distances between tree nodes and identification of internal MRCA nodes were performed using Bio.Phylo functions in Biopython (version 1.83).