The fresh leaves of A. odorata were collected in Sanya, Hainan, China. These specimens have been deposited in our lab (Seed Health Centre of China Agricultural University, Sanya Institute of China Agricultural University and Yunnan Agricultural University). Genomic DNA was extracted using the Tiangen Biotech DNA kit (Beijing). For library construction, we utilized the NEBNext^® library building kit with an insert size of 350 bp. The constructed DNA library was sequenced on the NovaSeq 6000 platform at Benagen (Wuhan, China). To ensure data quality, we applied Trimmomatic (Bolger et al., 2014) to remove low-quality sequences, including those with a quality value (Q) of less than or equal to 5, which accounted for more than 50% of the total bases, as well as sequences containing more than 10% “N” bases. Furthermore, the plant sample used for Illumina sequencing was also subjected to Oxford Nanopore sequencing based on PromethION devices. Purified DNA was prepared for long-read sequencing following the protocol outlined in the SQK-LSK109 genomic sequencing kit (ONT, Oxford, UK).

For long non-coding RNA (lncRNA) extraction, total RNA was isolated from fresh A. odorata leaves using a high-quality RNA extraction kit (TRIzol^® Reagent, Thermo Fisher Scientific, Waltham, MA, USA), following the manufacturer’s instructions. The extracted RNA was reverse-transcribed into cDNA using random primers, and rRNA was subsequently removed. The processed cDNA was fragmented and constructed into a library with an average length of 500 bp. The integrity and concentration of the RNA were assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The enriched lncRNA was then used to construct a cDNA library employing a protocol compatible with lncRNA sequencing. The library was subsequently sequenced using an Illumina HiSeq platform. Quality control measures were implemented to filter out low-quality sequences, and bioinformatics analyses were performed on the resulting data to identify and characterize lncRNAs.

For plastome assembly, we utilized GetOrganelle v1.7.4.1 with the following parameters: ‘-R 15 -k 21,45,65,85,105 -F embplant_pt’ to assemble the Illumina short-reads (Jin et al., 2020). GetOrganelle generated two complete plastome sequences, and we selected the one where the SSC region aligns in the same direction as Arabidopsis thaliana (NC_000932.1). Subsequently, we performed de novo assembly for long-reads. The long-reads were polished using Canu (Koren et al., 2017) and then assembled using PMAT assembler (Bi et al., 2024a) with the default parameters. BLASTn (Chen et al., 2015) was employed to identify the draft mitogenome from the assembled sequences. Six mitochondrial contigs were successfully identified. Considering the low accuracy of long-reads from Oxford Nanopore sequencing, we further utilize assembled mitochondrial sequences as the reference sequence. We establish an index using BWA (Li and Durbin, 2009). Subsequently, ‘bwa mem’ was utilized to obtain reads successfully mapped to the reference sequence. Indexing and mapping of the long-reads were performed using minimap2 (Li, 2018), with specific parameters ‘minimap2 -d’ for index creation and ‘minimap2 -ax map-ont -t 8 –secondary=no’ for mapping long-reads. Finally, we performed hybrid assembly using Unicycler (Wick et al., 2017) by combining Illumina short-reads and Nanopore long reads. The mapped Illumina short-reads were initially assembled using SPAdes (Bankevich et al., 2012), and then the Nanopore long-reads were employed to resolve repetitive sequence regions in the assembly, using minimap2 (Li, 2018). After multiple iterations and adjustments, we determined the optimal kmer value of 89. The resulting GFA format files generated by Unicycler were visualized using Bandage (Wick et al., 2015). Ultimately, Unicycler generated a complete circular genome. Notably, as these contigs were assembled based on Illumina short-reads, no additional polishing steps were necessary.

In our study, we employed PCR experiments to investigate the structure of A. odorata. Specifically, we designed eight specific primers to verify the accuracy of assembly generated by PMAT assembler. The primer design was conducted using the Primer designing tool on NCBI (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) with default parameters. The primer sequences used for PCR reactions are listed in Supplementary Table 1 . Subsequently, DNA was extracted, and the amplifications were performed using a Pro-Flex PCR system (Applied Biosystems, Waltham, MA, USA). The PCR reaction volume was 25 µL, comprising 2 µL of template DNA, 0.5 µL of forward primer, 0.5 µL of reverse primer, 12.5 µL of 2 × Taq PCR Master Mix, and 9.5 µL of ddH₂O. The amplification conditions consisted of an initial denaturation at 94 °C for 5 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 58 °C for 30 s, extension at 72 °C for 60 s, and a final extension step at 72 °C for 5 min. The PCR amplicons were visualized using 1% agarose gel electrophoresis. Subsequently, the single bright bands were excised and sent to Sangon Biotech (Shanghai, China) Co., Ltd. for Sanger sequencing.

The plastome of our A. odorata was annotated using CPGAVAS2 (Shi et al., 2019). The plastome of published A. odorata (NC_048994.1) was used as the reference genome. The annotation results were further verified using CPGView (Liu et al., 2023) to ensure accurate gene annotations. We utilized IPMGA (http://www.1kmpg.cn/ipmga/) to annotate the assembled mitogenome of A. odorata. We selected a database of mitochondrial genes of angiosperms on IPMGA. IPMGA generates annotated files in the standard GenBank format. The tRNA annotations were performed using tRNAscan-SE (Lowe and Eddy, 1997) while rRNA annotations were obtained through BLASTn (Chen et al., 2015). To ensure accuracy, manual edits were made to the annotations using Apollo (Lewis et al., 2002). Finally, the genome map was generated using OGDRAW (v1.3.1) (Alverson et al., 2010).

The long tandem repeats were detected by Tandem Repeats Finder (TRF, https://tandem.bu.edu/trf/trf.html) with the default parameters (Benson, 1999). The simple sequence repeats (SSRs) of the assembled mitogenome were identified using the online website MISA (https://webblast.ipk-gatersleben.de/misa/), the parameters of the minimum numbers of mono-, di-, tri-, tetra-, penta-, and hexanucleotides were set as 10, 5, 4, 3, 3, and 3, respectively. Additionally, forward, reverse, palindromic, and complementary repeat sequences were identified using REPuter (https://bibiserv.cebitec.uni-bielefeld.de/reputer/) with the following settings: hamming distance of three and minimal repeat size of 30 bp, and e-value is limited to less than 1e-5. The visualization of the repetitive elements was done using the Circos package (Zhang et al., 2013).

We employed PhyloSuite software (v1.2.2) (Zhang et al., 2020a) to parse the GenBank format file of the A. odorata mitogenome, extracting the protein-coding genes (PCGs). Subsequently, we conducted an analysis of the codon usage in mitochondrial PCGs using Mega 7.0 software (Kumar et al., 2016), which involved the calculation of Relative Synonymous Codon Usage (RSCU) values. An RSCU value of 1 signifies a neutral preference for codon utilization, whereas an RSCU value exceeding 1 indicates a relatively higher frequency of codon usage.

To identify the mitochondrial plastid DNAs (MTPTs), we compare the plastome and mitogenome DNAs of A. odorata by using BLASTn (Chen et al., 2015) program with the following parameters: -evalue 1e-5, -word_size 9, -gapopen 5, - gapextend 2, -reward 2, -penalty -3. The BLASTn results were visualized using Circos package (Zhang et al., 2013). The identified MTPTs were also annotated by using GeSeq.

We employed a two-step approach for predicting RNA editing sites. Initially, lncRNA-Seq reads were mapped to the coding sequences (CDS) of each protein-coding gene (PCG) using BWA software (Li and Durbin, 2009) with default parameters. Subsequently, we utilized REDItools (Picardi and Pesole, 2013) to predict RNA editing sites based on the mapping results. The prediction criteria were established as follows: coverage exceeding 30, frequency equal to or greater than 0.1, and p-value equal to or greater than 0.05. Afterward, the Illumina short-reads of DNA were aligned to the CDS of each PCG using BWA software with default parameters. Genomic SNPs were predicted using BCFtools (Li, 2011) based on the mapping outcomes, with thresholds set at coverage greater than 30 and frequency less than or equal to 0.1. These heterogeneous sites of natural variation need to be excluded from RNA editing sites. Finally, after excluding SNP sites, the remaining sites identified in the lncRNA-seq mapping will be considered as genuine RNA editing sites.

To confirm the accuracy of the predicted RNA editing sites, we further designed experiments to validate these four specific sites. The primers were designed on both sides of the editing sites ( Supplementary Table 2 ), and amplification was performed using genomic DNA (gDNA) and cDNA obtained from RNA reverse transcription using random primers as templates. The reaction conditions for amplification are described above. The amplified products were subsequently subjected to Sanger sequencing. Finally, by comparing the sequences of the products obtained from gDNA and cDNA, we determined the occurrence of RNA editing events.

For the collinear analysis with A. odorata, we selected five closely related species: Citrus unshiu (NC_057142.1), Citrus maxima (NC_057143.1), Citrus sinensis (NC_037463.1), Toona ciliata (NC_065060.1) and Toona sinensis (NC_065061.1). We identified collinear blocks based on sequence similarity using the BLASTn program with the following parameters: -evalue 1e-5, -word_size 9, -gapopen 5, - gapextend 2, -reward 2, -penalty -3. Only collinear blocks longer than 1 kb were retained for downstream analysis. To visualize the collinear relationships, we generated a multiple synteny plot using TBtools (Chen et al., 2023).

We retrieved a total of 31 mitogenomes, including two outgroups (Stylosanthes capitata and Glycine max), from the GenBank database. These mitogenomes were used to construct a phylogenetic tree with A. odorata. Firstly, PhyloSuite (v.1.2.2) (Zhang et al., 2020a) was employed to identify and extract orthologous protein-coding genes (PCGs) across the analyzed species. The nucleotide sequences corresponding to these PCGs were then aligned using MAFFT (v7.471) (Katoh and Standley, 2013). Subsequently, the aligned sequences were concatenated to generate the input for phylogenetic tree construction. The maximum likelihood (ML) method was implemented using IQ-TREE (version 2.1.4-beta) (Minh et al., 2020) with the parameters “–alrt 1000 -B 1000”. Using the Bayesian Information Criterion (BIC) for model selection, the results indicate that the best-fit model is GTR+F+R2. The bootstrap analysis was performed with 1,000 replicates. Finally, the resulting phylogenetic tree was visualized and edited using the online tool ITOL (Letunic and Bork, 2019).

The assembly is composed of six distinct nodes and eight edges (paths), visually depicted in Figure 1A . Each of these nodes signifies an assembled contig, demonstrating a region of overlap along the linkages. It’s worth noting that contig5 and contig6 exhibited distinct characteristics suggestive of potential repetitive sequences. These two repeat sequences each showcased four distinct paths (designated as p1–p4 and p5–p8, Figure 1B ). To confirm the presence of these paths within the A. odorata mitogenome, we conducted PCR experiments. The four primer pairs (F1 + R1, F2 + R1, F3 + R3, F3 + R2) were employed to validate the repeated sequence within contig5 and the paths, while the remaining pairs (F4 + R4, F5 + R5, F6 + R4, F7 + R7) were utilized to confirm the repetitive sequence within contig6 and the paths. The PCR products exhibited conformity with the anticipated outcome ( Supplementary Figure 1 ). The results of these PCR experiment not only validated the accuracy of the assembly and the eight paths, but also helped us to propose four possible genomic configurations. Configuration 1 presents a master circular structure, incorporating all six contigs ( Figure 1C ). In configuration 2, contig2, contig3, contig4, contig5, and contig6 collectively form a circular arrangement, and contig1 and contig6 form a smaller circular ( Figure 1D ). Similarly, configuration 3 showcases a circular arrangement encompassing all six contigs, while the contig3 was inverted compare to configuration 1 ( Figure 1E ). Lastly, configuration 4 is similar to configuration 2, with the inversion of contig3 ( Figure 1F ). Our PCR experiments show the possibility of multiple configurations, and here we use configuration 1, a master circle represent the complete mitogenome, for subsequent analysis.

The mitogenome maps of A. odorata was visually presents in Figure 2A . The total length of A. odorata mitogenome is 534,321 bp, and consists of 35 distinct protein-coding genes (PCGs) ( Table 1 ), including five ATP synthase genes (atp1, atp4, atp6, atp8, and atp9), four cytochrome c biogenesis genes (ccmB, ccmC, ccmFC, and ccmFN), nine NADH dehydrogenase genes (nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7, and nad9), three cytochrome c oxidase genes (cox1, cox2, and cox3), one transport membrane protein gene (mttB), one maturases gene (matR), and one cytochrome b gene (cob), four large subunits of ribosomal proteins (rpl2, rpl5, rpl10, and rpl16), five small subunits of ribosomal proteins (rps1, rps3, rps4, rps10 and rps12), as well as two succinate dehydrogenases (sdh3 and sdh4).

In the A. odorata mitogenome, a total of 22 tRNA genes have been annotated, with 19 being unique. Among these, 11 tRNA genes are mitochondrial native. Furthermore, our investigation has revealed 7 tRNA genes originating from the plastid: trnN-GUU, trnH-GUG, trnM-CAU, trnD-GUC, trnW-CCA, trnP-UGG, and trnI-CAU. Notably, our exploration has led us to the identification of a tRNA gene with bacterial origins, trnC-GCA, exhibiting a remarkable level of sequence homology with previously documented genes (Kitazaki et al., 2011). The remaining tRNA genes, devoid of sequence homology with known organelle tRNA genes, are unknown about their origin (Li, 2011; Rice et al., 2013). Furthermore, we have successfully pinpointed three distinct rRNA genes within the A. odorata mitogenome, namely rrn5, rrn18, and rrn26. The precise positions of each gene can be referenced in Supplementary Table 3 . Among the entirety of the genes that have been annotated, 10 PCGs encompass introns ( Supplementary Table 3 ). To elaborate, the genes ccmFC, cox2, cox1, trnP-CGG, and rps10 each feature one intron, whereas nad4 encompasses two introns. On the other hand, the genes nad1, nad2, nad5, and nad7 encompass four introns each.

Microsatellites, also known as simple sequence repeats (SSRs), typically consist of tandem sequences with a length of up to 6 base pairs in eukaryotic genomes. Within the mitogenome of A. odorata, a comprehensive total of 192 SSRs has been meticulously identified (as shown in Supplementary Table 4 ). Among this array of SSRs, tetrameric repeats stand out as the most prominent, encompassing 34.90% (67) of the overall count. This pattern is subsequently trailed by monomeric repeats (42), dimeric repeats (31), trimeric repeats (29), pentameric repeats (18), and hexametric repeats (5) ( Figure 2A ). And we have detected 29 long tandem repeat elements ( Supplementary Table 5 ).

In addition, we have meticulously identified a total of 333 pairs of dispersed repeats within the A. odorata mitogenome, each with lengths equal to or exceeding 30 base pairs. This collection encompasses 173 pairs of forward repeats, 146 pairs of palindromic repeats, and 7 pairs of reverse repeats and complementary repeats, as outlined in Supplementary Table 6 . Most of these repeat elements exhibit a length of less than 200 bp, and it’s evident that the number of dispersed repeats surpasses that of both SSRs and tandem repeats. The only two long dispersed repeats identified in the A. odorata mitogenome is contig5 and contig6 (10,113 bp palindromic repeat element and 6,621 bp repeat element). Cumulatively, the extent of these dispersed repeats encompasses 30,754 bp, constituting 5.76% of the entire A. odorata mitogenome. Furthermore, we provide a visual representation employing the Circos (v1.120) package (Krzywinski et al., 2009) to depict the dispersed repeats of the A. odorata mitogenome and counted the number of repetitive sequences in Figure 2B . These dispersed repeats are distributed in various regions of the mitogenome, effectively increasing the size of the genome.

We conducted an analysis of the codon usage within the PCGs. The comprehensive codon usage data for all PCGs is presented in Supplementary Table 7 . As shown in Figure 3A , revealing a discernible preference for specific codons among mitochondrial protein-coding genes. Notably, the RSCU values for the start codons AUG (Met) and UGG (Trp) both equate to 1. Furthermore, the RSCU values for the termination codons UGA (End), UAA (End), and UAG (End) are recorded as 1.20, 1.02, and 0.78, respectively. In terms of specific codons, GCU (Ala), UAA (End), CAU (His), CCU (Pro), and UAU (Tyr) emerge as the four most frequently employed codons within A. odorata. Conversely, GCG (Ala), UAG (End), CAC (His), and UAC (Tyr) are identified as the four least utilized codons. A visual representation in Figure 3A underscores the prevalence of arginine (Arg), leucine (Leu), and serine (Ser) codons, while methionine (Met) and tryptophan (Trp) codons exhibit relatively lower occurrence rates.

In our study, we conducted an annotation of the plastidial genome of A. odorata and performed a comprehensive comparison with its mitogenome. Employing the BLASTn program, we successfully identified a total of 12 instances of homologous sequences, we considered they might be potential MTPTs occurring between these two organelle genomes. These 12 MTPTs collectively span a length of 2,501 bp, contributing to 0.47% of the mitogenome’s total size ( Supplementary Table 8 ). Among these MTPTs, MTPT12 stands out as the longest, spanning 1,122 bp, while MTPT1 emerges as the briefest, encompassing a mere 29 bp. Subsequently, our efforts focused on annotating these MTPTs, revealing a revelation: every MTPT encompassed plastidial genes or gene fragments. As depicted in Figure 3B and Supplementary Table 8 , MTPT12 harbors a set of plastid genes, primarily associated with the photosystem II protein complex. It included psbJ, psbL, psbF, and psbE. Moreover, our analysis unearthed various gene fragments that resulted from the process of plastid migration. These fragments encompassed genes such as petG, ndhD, psbC, and atpH. It is plausible that these gene fragments underwent sequence loss during migration.

We have successfully discerned a total of 427 high-confidence C to U RNA editing sites across 32 mitochondrial protein-coding genes ( Supplementary Table 9 ). These editing sites are supported by our lncRNA data with an average depth of nearly 5100 times. The RNA editing sites for each gene are visually represented in Figure 3C . Within this set of mitochondrial genes, ccmB boasts the highest number of RNA editing sites at 45, closely followed by mttB with 41, positioning them as the foremost two genes in terms of RNA editing occurrences. Conversely, genes like atp1, atp6 and rpl5 possess the least number of editing sites, with only a solitary C to U edit detected for each. Based on our findings, we have identified C to U RNA editing events in three genes, which lead to the creation of premature stop codons. These genes are ccmFC, atp9 and rps10 (where CGA transitions to UGA). Notably, RNA editing plays a role in the formation of start codons as well, such as gene nad4L, have their start codons generated through RNA editing, converting ACG codon to AUG.

To evaluate the precision of this prediction, we employed PCR amplification and Sanger sequencing as a means of substantiating the occurrence of RNA editing. The four primer sequences can be found in Supplementary Table 2 . Among these, every RNA site underwent successful verification, namely nad4L-2, atp9-223, ccmFC-1324, and rps10-412 (as depicted in Figure 3D ). And all the validated sites encompassed C to U substitutions ( Figure 3E ), ccmFC-1324 is an exception, probably because of the lower frequency of editing, but the peak plot shows a hybrid peak of base U (uracil). Notably, the genomic DNA (gDNA) of rps10 has a longer PCR product because our primers span its introns, and the banding here also confirms that the DNA was sufficiently removed in the cDNA experimental group. Moreover, it’s noteworthy that each of these sites involved non-synonymous substitutions, as detailed in Supplementary Table 9 . These PCR experiment firmly attesting to the dependability of the anticipated RNA editing sites.

To delve into the rearrangements and conserved sequence blocks within the mitogenomes, we utilized the BLASTn program to pinpoint homologous collinear blocks. Illustrated in Figure 4A , each ribbon connecting two neighboring mitogenomes signifies a remarkably homologous collinear block or sequence. When comparing A. odorata and T. ciliata, we observed one large adjacent collinear block. Furthermore, we uncovered three collinear blocks exceeding 10 kb in length between these two mitogenomes. However, in the case of comparing A. odorata with C. sinensis, no extensive collinear blocks were detected, with the longest identified block measuring only 7.5 kb. Overall, while longer collinear blocks tend to be present between closely related species, the mitogenomes demonstrated limited collinearity, with several regions lacking homology. These findings underscore the prevalence of extensive genomic rearrangements between A. odorata and its related mitogenomes, suggesting that the genomic structure of the mitogenomes is not conserved.

In addition, we carried out a phylogenetic analysis utilizing 31 mitogenomes of related species, with T. capitata and G. max serving as outgroups for reference. The detailed list of species and their corresponding GenBank accessions utilized for this analysis is available in Supplementary Table 10 . By aligning and concatenating the shared protein-coding genes (PCGs), we created the matrix data for analysis. The outcome of our phylogenetic investigation yielded a maximum likelihood (ML) tree that exhibits robust support along the primary basal branches (refer to Figure 4B ). In terms of evolutionary relationships, A. odorata exhibits its closest affinity with T. ciliata, both species falling under the Sapindales order. While the overall structure of the phylogenetic tree broadly aligns with the APG IV system (The Angiosperm Phylogeny et al., 2016) at the order level, it’s noteworthy that within the Sapindales order, a couple of nodes lack bootstrap support. This suggests that deriving phylogenetic inferences solely from mitochondrial PCGs might have limitations in accurately resolving lower taxonomic categories.