PacBio Sequel-II raw subread data generated from the whole genome of E. camaldulensis was downloaded from NCBI and was used for the mitogenome assembly (SRX11929917). High Fidelity Circular Consensus reads (HiFi CCS) were locally called from the downloaded subread data using SMRTLink ver. 11.1. The first de novo assembly was computed using hifiasm ver. 0.16.1 (Cheng et al., 2021) against the whole PacBio HiFi dataset. Hereafter we refer to contigs computed by hifiasm as those computed with the whole-genome sequencing dataset, including nuclear genome assembly (Driguez et al., 2021). We searched mitogenome contigs with homology to the mitogenome of E. grandis (NCBI reference sequence: NC_040010.1). Because the collected contigs show different forms of the sequences, rearrangements mediated by repeat sequences were considered to explain multiple contigs (Sloan, 2013). Therefore, a repeat-aware graph assembler was applied to suppress this issue (Fischer et al., 2022). To accelerate computational time, mitochondria-derived reads were extracted from the original dataset using the E. grandis genome with minimum query coverage to the reference at 70% using minimap2 ver.2.26 (Li, 2018). We assembled the E. camaldulensis mitogenome using Flye ver. 2.8.3 with default parameters (Kolmogorov et al., 2019). The assembled mitogenome contig was annotated by using MITOFY ver. March 22, 2012 (Alverson et al., 2010) and tRNAscan-SE ver. 2 (Chan et al., 2021) with manual corrections to confirm the annotation. The genome circle map was drawn by using ORGDRAW ver. 1.3.1 (Greiner et al., 2019). The mitochondrial reads were reextracted from the entire dataset by mapping against reference sequences containing our mitogenome assembly, the nuclear genome (GCA_019915185.1), and the plastid genome (CM024681.1) for further analysis.

The direct repeat (DR) and inverted repeat (IR) contents of the E. camaldulensis mitogenome were analyzed using ROUSfinder ver. 2.0 and blastn ver. 2.14.1 (Altschul et al., 1997; Wynn and Christensen, 2019). Validation of the repeat-mediated genomic rearrangement was performed by mapping the mitochondria-derived reads against all possible sequences with the target repeat in the middle. Specifically, for each repeat pair detected by ROUSfinder, 3,000-bp flanking sequences upstream and downstream of the repeat sequence were extracted, and four sequences representing possible isomers mediated by homologous recombination were generated with the extracted flanking and the target repeat sequences. These generated sequences were combined as a reference set. The mitochondria-derived reads were mapped to the reference sequence set using minimap2 ver. 2.26, and supporting reads were counted for each repeat pair if a read is mapped to a sequence with at least 1,800-bp overlap to flanking sequences in both ends as a primary alignment. A few repeats are nested in another larger repeat sequence; however, the 1,800-bp threshold is long enough to detect unique sequences beyond the repeat sequence covering the target nested repeat (the largest covering repeat sequence has 1,467 bp in length). The assembly used for the repeat analysis was named Master Circle 1 form (MC1). It should be noted that because the rates of rearrangement mediated by the long (>1,000 bp) repeat pairs were significant, we designated the rearranged MC1 sequence mediated by the longest repeat (repeat-1) as Master Circle 2 form (MC2). Additionally, the sequences of MC1 and MC2 rearranged by the second longest repeat (repeat-2) were termed MC1² and MC2², respectively.

The mitogenomes of E. grandis and E. camaldulensis were compared using LAST ver. 1454 (Kiełbasa et al., 2011) with default parameters. The MC2 of the E. camaldulensis mitogenome was selected as a representative sequence due to the least number of variations compared to the E. grandis mitogenome. Detected significant hits were visualized and summarized using the maf-convert command. Large rearrangements were visualized by using SyRI (Goel et al., 2019). An E-value of 1e-20 was used as the threshold.

The taxonomic placement of E. camaldulensis was established by analyzing the complete mitochondrial genomes of all plant species within the Myrtaceae family available in the NCBI database. The complete genomes of E. grandis, Syzygium samarangense, and Rhodomyrtus tomentosa (accession numbers: NC_040010.1, NC_079700.1, and NC_071968.1, respectively) were downloaded. The complete mitochondrial genome of Begonia coptidifolia (accession number: LC706752.1) was used as an outgroup for rooting purposes. A multiple sequence alignment of five mitochondrial genomes, including E. camaldulensis, was generated using progressiveMauve with default parameters (Darling et al., 2010). Phylogenetic relationships were reconstructed using maximum likelihood with the general time reversible model. Initial trees for the heuristic search were obtained automatically by applying neighbor-joining and BioNJ algorithms to a matrix of pairwise distances estimated using the maximum composite likelihood approach and then selecting the topology with a superior log likelihood value. To model differences in the rates of evolution between sites, a discrete Gamma distribution was used. Positions containing missing/ambiguous data were excluded. An estimate of support for each branch was determined with the bootstrap test using 100 pseudoreplications. Phylogenetic analyses were conducted in MEGA11 (Tamura et al., 2021).

Fresh E. camaldulensis leaves were collected from the research greenhouse at King Abdullah University of Science and Technology, Thuwal, Saudi Arabia. HMW DNA was extracted by following a modified Qiagen protocol (dx.doi.org/10.17504/protocols.io.bafmibk6).

This study utilized the GridION sequencer for Oxford Nanopore sequencing of E. camaldulensis. In preparation for library creation, short fragments were removed from the extracted HMW DNA using either BluePippin with a 30-kb cutoff or the Short Read Eliminator XL kit (SRE XL, Circulomics, Baltimore, MD, USA) as per the manufacturer’s instructions. The size-selected DNA underwent a cleanup process with AMPure XP beads, followed by library processing using the genomic DNA by ligation protocol (SQK-LSK109, Oxford Nanopore, Oxford, UK). Specifically, DNA samples of E. camaldulensis (1.3 μg) were repaired and end-prepped before bead cleanup and adapter ligation. The ligated product underwent bead cleaning with long fragment buffer (LFB) and was subsequently eluted. The prepared library for E. camaldulensis (8 fmol) was loaded onto the GridION platform (FLO-MIN106D). This study extracted mitochondrial reads from the ONT dataset using the same criterion applied for the HiFi dataset. The duplex-like reads (i.e., reads mappable to the same region on the reference in two strands) were removed because they cannot be mapped lineally.

The mitochondrial reads in both HiFi and ONT datasets were mapped against MC1 sequence. Only the reads that had 95% of their bases aligned linearly with the MC1 sequence were retained for coverage computation. Reads that were partially mapped, i.e., those with less than 95% of their bases were aligned linearly, were mapped against an alternative sequence, MC2². This was done if a portion of the reads fully and linearly aligned with the alternative form but not with the reference sequence. To filter reads, CoverM ver. 0.7.0 was utilized (https://github.com/wwood/CoverM).

The assembled E. camaldulensis mitogenome is a single circular sequence of 463,134 bp, with an overall GC content of 45.11% ( Figure 1 ). The mitogenome has a quadripartite structure that is characterized by a pair of long (8,986 bp) IR sequences ( Supplementary Figure S1 ). Our findings suggest that a circular form, rather than the hypothetical linear structure previously discussed in the case of E. grandis (Pinard et al., 2019), is dominant in the mitogenome of E. camaldulensis. In our E. camaldulensis assembly, we did not find any results that support the linear structure. However, it should be noted that this does not exclude the possibility that the mitogenome contains a head-to-tail concatemer structure in vivo (Sloan, 2013).

Annotations were made for 38 conventional protein-coding genes within the single circle model. In addition, 19 tRNA genes and five rRNA genes were observed ( Table 1 ). The E. camaldulensis mitogenome seems to have lost rps2 and rps1 genes (Adams et al., 2002). Although fragmented homologous sequences of rps19 were detected in the mitogenome, stop codons were also observed in the middle of the gene. Two copies of 5S rRNA were found, which are also present in the mitogenome of E. grandis. Two copies of 18S rRNA were also annotated on the E. camaldulensis mitogenome, which was previously described as a single copy gene in E. grandis mitogenome (Pinard et al., 2019). Overall, it appears that the gene content is well conserved between the two Eucalyptus mitogenomes.

To date, there are only three complete mitochondrial sequences from the Myrtaceae family being publicly available. For the sake of completeness, we performed a phylogenetic analysis to ascertain the phylogenetic placement of the E. camaldulensis mitogenome. The topology of the phylogenetic tree suggests a strong phylogenetic connection between the two eucalypt species, E. grandis and E. camaldulensis ( Figure 2A ). The tree topology aligns with the phylogenetic tree inferred using mitogenomes by Lu and Li (2024).

The co-transcription and physical proximity of 18S rRNA and 5S rRNA in plant mitogenomes have been reported (Chen et al., 2017). Given the topology of the phylogeny, the varying copy numbers of these genes suggest multiple hypotheses for the evolution of these two genes on the mitogenome in the Myrtaceae family ( Figure 2B ). The most parsimonious hypothesis is that a single duplication of 5S rRNA occurred in a common ancestor of the two eucalypts, followed by a duplication of 18S rRNA in the lineage of E. camaldulensis. This hypothesis assumes that there was a duplication event of 18S rRNA into a proximal physical location of 5S rRNA. Alternatively, another hypothesis suggests that there were two duplication events of the two genes in a common ancestor, with a subsequent loss of 18S rRNA that happened in the lineage of E. grandis. Estimating the most likely evolution requires additional information regarding the state of the last common ancestor and more complete mitogenomes.

The number of detected dispersed repeat sequences was 63, and the repeat size ranged from 50 bp to 9 kb ( Supplementary Table S1 ). In the long repeat category (>1,000 bp), we identified pairs of IRs and DRs with lengths of 8,986 and 1,467 bp, respectively. As it has been argued that rearrangement mediated by long repeats is frequent in the plant mitogenomes (Siculella et al., 2001; Sloan, 2013), other possible isomeric forms were investigated by read mapping (see “Materials and methods”). Rearrangements by long repeat sequences are more common than those mediated by middle or short repeat sequences ( Table 2 ). In particular, isomeric forms mediated by the pair of the largest IR seem to exist at an equal frequency, similar to other plant mitochondria (Sloan et al., 2010) and chloroplast genomes (Wang and Lanfear, 2019). Rearrangement mediated by the second largest DR pair is also common. The data are consistent with a multipartite mitochondrial genome as a two-circle model ( Figure 3A ). As an attempt to validate these rearrangements mediated by dispersed repeats, we analyzed our assembly using the Oxford Nanopore sequencing platform as an alternative and independent experiment (“Materials and methods”). Using the two independent sequencing technologies, we were able to validate the rearrangements mediated by repeat 1 and repeat 2 ( Figure 3B ). There are also other rearrangements caused by shorter repeats (100 bp to 1 kb), but their frequency is significantly lower ( Table 2 ; Supplementary Table S2 ). In fact, 95% of mitochondrial reads mapped against either MC1 or MC2² fully and linearly in both the HiFi and ONT sequencing datasets, and the rest of the mitochondrial reads are mainly related to the minor rearrangements ( Supplementary Table S3 ). The HiFi sequencing data and further validation by ONT altogether support four major isomeric forms, including the single circle model ( Figure 3 ).

Since the contigs computed by hifiasm only partially matched the major isomeric sequences, we investigated the multiple contigs generated by hifiasm to understand the technical challenges for plant mitogenome assembly. The size of repeat-5 is 188 bp in length, and it is a part of another larger repeat-2 that has a length of 1,467 bp (188 bp of 1,173-th to 1,360-th positions of repeat-2). Repeat-2 is related to one of the major rearrangements ( Figure 3 ), and repeat-5 is detected by a 188-bp long sequence far away from the repeat-2 sequences ( Supplementary Table S1 ). Rearrangements explained by those small and nested repeat pairs appear to be important despite the low frequency ( Table 2 ). Those minor isomeric sequences aligned well to the unique sequences found in the contigs computed by hifiasm, a haplotype-aware assembler ( Figure 4 ; Supplementary Figure S2 ). A haplotype-aware assembler such as hifiasm would have concatenated major and minor isomeric forms together as single haplotypes ( Figure 4B ), and this seems to be the cause of multiple contigs.

A comparative analysis was performed on the mitogenomes of E. camaldulensis and E. grandis to understand how synteny is preserved between these two mitogenomes within the same genus. Overall synteny was found to remain intact except for a large inversion observed between positions 18,587 and 224,548 in the E. grandis mitogenome (corresponding to positions 19,028 and 212,830 in E. camaldulensis) ( Figure 5A ). Despite the fact that synteny is still maintained to some degree, it does not explain a major difference: the hypothetical linear and circular form for them. Due to the absence of other reported mitogenomes in the Eucalyptus genus, phylogenetic reconstruction of the ancestral form remains infeasible. However, intriguing fragments with homology to the 5′ terminal of the E. camaldulensis genome were detected at both the 5′ and 3′ terminals of the E. grandis mitogenome ( Figure 5B ). As indicated earlier, while the E. camaldulensis mitogenome encodes two copies of 18S rRNA ( Table 1 ), only one copy is retained in the E. grandis mitogenome ( Figure 2B ). Given that this gene is encoded within repeat-1 of the E. camaldulensis mitogenome ( Figure 5B ), we can envisage that this gene could serve as a crucial marker for historical tracking without necessitating additional mitogenomes. A sequence comparison using LAST unveiled mutated and fragmented sequences of 18S rRNA in the E. grandis mitogenome ( Figure 5B ). Notably, a fragment resembling a second half of the gene was located at the very 5′ end of the E. grandis mitogenome, and a fragment homologous to the first half of the gene was also recognized at the opposite end (3′ end) ( Figure 5B ). It is worth noting that the E. grandis mitogenome has unique and non-homologous flanking sequences at both terminals, which is consistent with the genome’s proposed linear structure.

We identified four repeat pairs—DR-1 (1,148 bp), DR-2 (366 bp), DR-3 (4,210 bp), and DR-4 (206 bp)—in E. grandis mitogenome (Pinard et al., 2019) as homologous to repeat-1, the largest IR in E. camaldulensis ( Figure 5 ). The combined length of these four repeat pairs is still shorter than that of repeat-1 (5,930 vs. 8,986 bp), suggesting that a split event may have occurred. On the opposite side of the E. grandis genome, a similarity to repeat-1, identified as fragmented 18S rRNA, is also observed ( Figure 5B ). Assuming that the homologous fragments at both the 5′ and 3′ ends of this mitogenome originated from a single repeat sequence (i.e., repeat-1a), it is possible that either fragment could retain similarity to the homolog of its counterpart, termed repeat-1b. Pinard et al. described three structural variations, termed BP1, BP2, and BP3, with a high degree of confidence (Pinard et al., 2019). Of these variations, BP2 was associated with a rearrangement event triggered by a pair of DRs observed in the E. grandis mitogenome. However, no such associations could be made for the other two variations, BP1 and BP3. Notably, the breakpoints of BP3 coincided with two homologous regions 2 kb distant from the start of repeat-1 ( Figure 6 ). Furthermore, one breakpoint of BP1 overlapped with a region at the start of repeat-1b, and the other breakpoint coincided with another region homologous to the end of the E. camaldulensis mitogenome. Due to its continuous nature (represented as circular form), this end region seamlessly transitions into the start of repeat-1a in the E. camaldulensis mitogenome. These structural variants together show a strong correlation (within 250 bp) with respect to the E. camaldulensis coordinates ( Figure 6 ).