Young tender leaves of I. indigotica were sourced directly from Xianglian Horticulture, a verified commercial nursery located in Weitang Town, Xiangcheng District, Suzhou, China (31.2250° N, 120.6420° E) ( Figure 1A ). The identity of the entire plant was initially verified in-house by comparing morphological characteristics to the reference specimen image IMC0058462 from the National Plant Specimen Resource Center of China Digital Herbarium (NPSRCCDH). To reinforce the accuracy of this species determination, two independent experts from Jiangsu University—both with extensive experience in plant taxonomy—were subsequently consulted and confirmed the identification. Young leaves were freshly collected directly from potted I. indigotica plants and immediately rinsed with DEPC-treated water. The cleaned tissue was flash-frozen in liquid nitrogen and stored at −80 °C to preserve DNA integrity. High-molecular-weight genomic DNA was extracted using a modified CTAB protocol optimized for plant materials. DNA quality was assessed by 0.75% agarose gel electrophoresis, while purity and concentration were determined using both UV spectrophotometry (NanoDrop One, Thermo Fisher Scientific) and fluorometry (Qubit 3.0, Life Technologies). Only high-quality DNA was used for library construction with the SMRTbell Express Template Preparation Kit 2.0 (Pacific Biosciences), followed by long-read sequencing on the PacBio Sequel II system for accurate organelle genome assembly.

The mitogenome of I. indigotica was assembled using a multi-stage pipeline optimized for PacBio HiFi reads. To identify mitochondrial-derived sequences, HiFi reads were first aligned against three published Brassicaceae mitochondrial reference genomes: Brassica rapa (NC_049892.1) (Ren et al., 2023), R. sativus (NC_018551.1) (Chang et al., 2013), and Sinapis arvensis (KM851044.1) (Sang et al., 2020). Prior to alignment, BLAST databases were constructed for each reference using makeblastdb from BLAST+ (v2.13.0+), with default parameters and nucleotide indexing (Camacho et al., 2009). Read alignment was then conducted using blastn with the following parameters: -evalue 1e-5 -outfmt 6 -max_hsps 10 -word_size 7 -task blastn-short. The resulting alignments were filtered to retain only uniquely matched reads, and read IDs from all three comparisons were merged to generate a non-redundant mitochondrial candidate set. These filtered reads were extracted using seqkit grep and subsequently assembled de novo with PMAT (Plant Mitochondrial Assembler Tool), employing the autoMito module with -st hifi mode, an estimated genome size of 300 Mb, and 24 computational threads. To improve the accuracy of the long-read data, PMAT, through its “correct_sequences.py” script, calls the Canu program to perform error correction, utilizing overlap-based assembly algorithms and error correction techniques. This procedure yielded a complete, circularized mitogenome of I. indigotica suitable for downstream analysis. To evaluate the assembly structure, the graphical layout in GFA format was inspected using Bandage (v0.8.1), allowing visual confirmation of the circular topology and connection paths among contigs (Wick et al., 2015).

After genome assembly, structural gene annotation was carried out using the Plant Mitochondrial Genome Annotator (PMGA, v1.5.3), a specialized tool designed for organelle genome annotation workflows (Li et al., 2025). The software was deployed within a reproducible computational environment managed by mamba, and executed inside a containerized system based on singularity (v3.7.2), ensuring compatibility and scalability across platforms. To improve annotation fidelity, the internal reference dataset (Dataset1) bundled with PMGA was utilized, offering curated models of mitochondrial genes. The initial automated annotation output was carefully reviewed and refined using Apollo (v2.5.0) (Lee et al., 2013), with manual correction of ambiguous gene boundaries and confirmation of structural features. The resulting high-confidence GFF3 annotation file, after manual curation, was submitted along with the corresponding FASTA genome sequence to the GenBase database (https://ngdc.cncb.ac.cn/genbase/?lang=en) for public access and future reference (Bu et al., 2024).

Repetitive sequences in eukaryotic genomes are broadly classified into simple sequence repeats (SSRs), tandem repeats, and dispersed repeats based on their structural features and genomic organization (Han et al., 2022; Ala et al., 2023; Li et al., 2023; Zhang et al., 2023; Wang et al., 2024). In this study, SSRs in the mitogenome of I. indigotica were identified using the misa.pl script (v2.1) with predefined thresholds: a minimum of 10 repeat units for mononucleotides, 5 for dinucleotides, and 4, 3, 3, and 3 for tri-, tetra-, penta-, and hexanucleotides, respectively (Thiel et al., 2003). A maximum distance of 1,000 bp was allowed between neighboring SSRs to detect potential compound repeats. Tandem repeats were identified using Tandem Repeats Finder (TRF, v4.09) under the parameter setting “2 7 7 80 10 50 500 -f -d -m”, which is optimized for the detection of tandem motifs in organelle genomes (Benson, 2023; Zhang et al., 2024). Dispersed repeats, defined as non-contiguous homologous sequences, were detected using the REPuter tool with the following parameters: minimum repeat length of 30 bp, maximum of 500 repeats per genome, and a Hamming distance threshold of 3 to accommodate minor mismatches (Kurtz et al., 2023). To gain insights into the distribution and organization of repeat elements across the mitogenome, all detected SSRs, tandem repeats, and dispersed repeats were mapped using Circos software (v0.69-8) (Krzywinski et al., 2023).

To assess whether dispersed repeat pairs mediate homologous recombination, two structural models were constructed for each candidate pair: (i) a reference configuration representing the native genomic arrangement, and (ii) a recombinant configuration simulating the outcome of recombination (Gualberto et al., 2014; Gualberto and Newton, 2017). For each model, 1000 bp flanking sequences surrounding the repeat pair were extracted. Nanopore long reads were then mapped to both models, and support for the recombinant structure was considered indicative of active recombination at that locus.

To experimentally validate selected recombination sites, PCR amplification was performed using primers designed to span the flanking regions of the repeat pair. Specifically, primers were designed in the 1000 bp regions upstream and downstream of each repeat: a forward primer (F) in the upstream region and a reverse primer (R) in the downstream region. For dispersed repeat pairs that are located far apart in the mitogenome, two sets of primers were used: F1/R1 and F2/R2 for the two respective repeat loci. If direct (forward) recombination occurs, it can be detected by amplification of bands using the F1/R2 and F2/R1 primer combinations. Conversely, if inverted recombination occurs, it may be evidenced by amplification with F1/F1 or R1/R2 combinations. These diagnostic PCR patterns provide direct molecular evidence for the type and presence of recombination events mediated by repeat sequences [refer to Figure 2 in (Liu et al., 2025)].

Each 25 μl PCR reaction contained 1 μl genomic DNA, 1 μl each of 8 μM forward and reverse primers, 13 μl of 2× Taq PCR Master Mix, and 10 μl nuclease-free water. Thermocycling conditions were as follows: 96 °C for 3 min initial denaturation; 32 cycles of 94 °C for 30 s, 62 °C for 30 s, and 72 °C for 1 min; and a final extension at 72 °C for 10 min. PCR products were analyzed by agarose gel electrophoresis, and positive bands were sequenced using the Sanger method to confirm recombination events and evaluate repeat-associated structural variation.

During the structural correction phase of mitogenome assembly, the complete cpgenome of I. indigotica was incidentally resolved as an independent circular contig within the graphical assembly output. This structure was identified and separated using Bandage (v0.8.1) through examination of the GFA-format assembly graph generated from long-read sequencing data ( Figures 1B,C ). The circular chloroplast sequence was subsequently exported in FASTA format for further characterization.

Annotation of the cpgenome was conducted using the GeSeq web server, with I. indigotica Fortune employed as the reference genome to guide gene prediction and structural annotation (GenBank: KT939360) (Yang and Wang, 2017). The annotated cpgenome was visualized using OGDRAW (v1.3.1) (Greiner et al., 2019), resulting in a high-resolution circular genome map. The generated GFF3 file was subsequently curated manually to correct annotation discrepancies and feature boundaries. Both the finalized GFF3 and the corresponding FASTA file were submitted to the GenBase database for long-term archival and community access.

Species-level validation of the cpgenome was performed through comparative analysis with two previously published plastomes: I. indigotica (GenBank: OP620952) and I. tinctoria (GenBank: NC_028415), confirming its taxonomic identity.

To investigate potential sequence transfer events between genomic compartments in I. indigotica, a comprehensive analysis of mitochondrial-plastid DNA transfer (MTPT) was conducted. The analysis employed BLASTn (v2.13.0) with the plastid genome used as the query and the mitogenome as the target. A stringent e-value cutoff of 1e-6 was set to ensure the detection of only high-confidence homologous regions (Chen et al., 2015). Regions of plastid origin identified within the mitogenome were visualized using TBtools (v2.010). The corresponding sequences were further annotated via GeSeq to assign functional roles and validate their organellar origin (Chen et al., 2020).

In parallel, the nuclear genome of I. indigotica (GenBank accession: VHIU00000000), which comprises seven chromosomes in its haploid complement, was also analyzed to assess possible organelle-to-nucleus DNA integration. This tri-compartmental comparative approach enabled the identification and characterization of both plastid-to-mitochondria and organelle-to-nucleus horizontal sequence transfers, offering insights into the dynamic genome interactions and evolutionary history of the species.

RNA editing is a crucial post-transcriptional process in plant mitochondria, predominantly characterized by cytidine-to-uridine (C-to-U) conversions. These events can restore conserved amino acids, generate functional start or stop codons, and impact mitochondrial protein function and stability (Lukes et al., 2021). To comprehensively identify RNA editing events in I. indigotica, we employed two complementary strategies: (1) RNA-seq read mapping and variant calling, and (2) comparative prediction using conserved reference genomes.

For the mapping-based approach, publicly available RNA-seq data (SRA: SRR9329298, BioProject: PRJNA549758) were aligned to a custom-assembled mitogenome of I. indigotica (genbase: C_AA108663) (Kang et al., 2020). The processing pipeline consisted of fastp (v0.24.0) for adapter removal and quality filtering, fastQC (v0.12.1) for quality assessment, HISAT2 (v2.2.1) for splice-aware alignment, and SAMtools (v1.6) for BAM file handling (Li et al., 2009; Kim et al., 2015; Chen et al., 2018; Wingett and Andrews, 2018). Variant calling was performed using BCFtools (v1.21), and RNA editing candidates were filtered using custom scripts developed in-house (Narasimhan et al., 2016; Danecek and McCarthy, 2017; Genovese et al., 2024). In parallel, a reference-based prediction strategy was implemented using the PREPACT3 platform (http://www.prepact.de/prepact-main.php, accessed on 26 Jun 2025; version 3.12.0/2.2.26+) (Lenz et al., 2010). To maintain high specificity, only two mitochondrial reference genomes were used: Brassica napus (NC_008285.1) and Citrullus lanatus (NC_014043.1). A stringent BLASTX e-value threshold of 1e−3 was applied to identify conserved editing sites with high confidence, following established methodologies (Lenz and Knoop, 2013).

Codon usage optimization—defined as the preferential employment of synonymous codons to improve translational fidelity and effectiveness—is a recognized phenomenon in both unicellular and multicellular organisms (Presnyak et al., 2015; Sen et al., 2020). This codon usage bias arises from a combination of selective pressures, species-specific translational demands, and random genetic events (Hanson and Coller, 2018). In this study, mitochondrial protein-coding proteins (PCGs) of Eleutherococcus senticosus were subjected to codon bias evaluation. The analysis was conducted using TBtools software (v2.010), and relative synonymous codon usage (RSCU) values were determined via CodonW (v1.4.4) (Charif et al., 2005).

Assessing the evolutionary patterns of PCGs is essential for understanding molecular adaptation and divergence among closely related species. In this study, the rates of nonsynonymous (Ka) and synonymous (Ks) substitutions were estimated to evaluate the selective pressures acting on mitochondrial PCGs of I. indigotica in comparison with 15 other representative species from the order Brassicales. The Ka/Ks ratio serves as an indicator of selection: values >1 suggest positive selection, = 1 imply neutral evolution, and <1 indicate purifying selection.

Homologous gene pairs between I. indigotica and the other Brassicales species were identified using BLASTN (v2.10.1). Shared mitochondrial PCGs were aligned using MAFFT (v7.313) in auto mode to ensure accurate multiple sequence alignment (Rozewicki et al., 2019). Ka and Ks values were computed with the Maximum Likelihood (MLWL) method implemented in Ka/Ks Calculator v2.0 (Zhang et al., 2006; Wang et al., 2010; Zhang, 2022). The distribution of Ka/Ks ratios across genes was visualized using boxplots generated with the R package ggplot2. In parallel, to assess genetic variability at the nucleotide level, aligned homologous gene sequences were analyzed using DnaSP v5, and nucleotide diversity (Pi) values were calculated for each PCG (Rozas et al., 2003; Librado and Rozas, 2009; Rozas et al., 2017). This analysis provides insights into the extent of sequence polymorphism within mitochondrial genes, further contributing to our understanding of evolutionary constraints and divergence patterns across Brassicales mitogenomes.

To elucidate the phylogenetic position of I. indigotica within angiosperms, a comparative mitogenomic analysis was performed using mitogenome sequences from 26 plant species representing four taxonomic orders. One species from the order Solanales was included as the outgroup to root the phylogenetic tree ( Supplementary Table S1 ). Among the selected species, 16 belonged to the family Brassicaceae, including I. indigotica itself. These cruciferous taxa were: Arabidopsis thaliana, Arabidopsis lyrata, Boechera stricta, Brassica napus, B. rapa, Brassica oleracea, Brassica juncea, Brassica carinata, Capsella rubella, Crucihimalaya lasiocarpa, Descurainia sophia, I. tinctoria, Lepidium apetalum, Lepidium sativum, and R. sativus. To provide broader phylogenetic context, 10 additional non-Brassicaceae species were included: Apium graveolens (Apiaceae), Panax quinquefolius (Araliaceae), Phaseolus vulgaris, Aeschynomene indica, Senna tora, and Vigna angularis (all Fabaceae), Saussurea costus and Taraxacum mongolicum (Asteraceae), and the outgroup Solanum aethiopicum (Solanaceae). This diverse taxon set enabled robust resolution of both intra- and inter-family evolutionary relationships.

All genomic datasets were retrieved from GenBank and subsequently processed using PhyloSuite (v1.2.3), which facilitated the extraction and reformatting of mitochondrial coding sequences (Zhang et al., 2020; Xiang et al., 2023). Multiple sequence alignment (MSA) of conserved gene regions was performed using MAFFT (v7.313) to ensure high-accuracy alignment across all taxa involved (Rozewicki et al., 2019). To determine the most appropriate nucleotide substitution models for downstream phylogenetic inference, PartitionFinder2 was employed. Maximum likelihood phylogenetic reconstruction was carried out using IQ-TREE2 (v2.1.4) with ultrafast bootstrap support (Minh et al., 2020). The selection of the best nucleotide substitution models was based on the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), which evaluated various models to identify the one that best fit the data. The final result indicated that the GTR+F+I+G4 model was the most appropriate, as it provided the lowest log-likelihood value and the best fit according to both AIC and BIC criteria. The final tree visualization and annotation were achieved through iTOL (v6), offering customizable graphics for enhanced biological interpretation (Letunic and Bork, 2021).

To investigate mitogenome collinearity and structural conservation in I. indigotica, three additional mitogenomes—Apium graveolens, Taraxacum mongolicum, and Vigna angularis—were selected for comparative analysis. These species, representing different families within the eudicots, provided a broad phylogenetic context for assessing syntenic relationships. Pairwise comparisons were conducted using BLASTn, with an e-value cutoff of 1e-6 to ensure specificity. Homologous regions exceeding 500 bp in length were designated as conserved collinear blocks. The collinearity relationships among I. indigotica, I. tinctoria, and the three selected eudicot species were visualized using NGenomeSyn (v1.41), enabling detailed inspection of genome rearrangements and conserved structural segments across lineages (He et al., 2023).

To assess large-scale sequence conservation and structural variation among closely related cruciferous species, comparative alignment of mitogenomes from I. indigotica, I. tinctoria, and R. sativus was conducted. Whole-genome alignments were generated using long-read-assembled mitogenomes and visualized with Mauve (v2.4.0) alignment software (https://sourceforge.net/projects/mauve). The Mauve viewer enabled the detection of locally collinear blocks (LCBs) and facilitated inspection of genome rearrangements, inversions, and insertions/deletions across the three Brassicaceae species.

A total of 368,520 HiFi long reads were generated using the PacBio Revio platform, producing approximately 6.48 Gb of high-quality sequencing data. The mean read length was 17,596 bp, with a median length of 16,546 bp and a read length N50 of 17,479 bp, indicating a highly consistent and uniform long-read distribution across the dataset ( Supplementary Figure S1A ). The median read quality reached Q33, demonstrating the exceptional base-level accuracy achieved through PacBio’s circular consensus sequencing (CCS) technology ( Supplementary Figure S1B ). The sequencing depth for the mitochondrial genome was approximately 24.9x, and for the chloroplast genome, the sequencing depth was approximately 42.2x, ensuring comprehensive coverage of both genomes. To ensure data transparency and facilitate future reuse, the raw HiFi sequencing data have been deposited in the Genome Sequence Archive (GSA) at the National Genomics Data Center (NGDC) under accession number CRA026873 (https://ngdc.cncb.ac.cn/gsa) (Members and Partners, 2025). Furthermore, the annotated mitogenome assembly has been submitted to GenBase, also hosted by NGDC/CNCB, and is accessible under accession number C_AA108663.1 (https://ngdc.cncb.ac.cn/genbase) (Bu et al., 2024). The final assembly comprises a complete circular mitogenome of 260, 864 base pairs in length ( Figure 1D ).

The mitogenome of I. indigotica, as annotated using the PMGA pipeline, comprises 65 unique genes, including 23 core PCGs, 9 variable genes, 3 ribosomal RNA (rRNA) genes, and 18 distinct transfer RNA (tRNA) species, with copy number variation observed in several tRNAs ( Table 1 ). The core PCGs are involved in essential mitochondrial functions and include five ATP synthase subunits (atp1, atp4, atp6, atp8, atp9), nine NADH dehydrogenase genes (nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7, nad9), one ubiquinol cytochrome c reductase gene (cob), and six genes related to cytochrome c biogenesis, with duplicated copies of ccmB and ccmFN (ccmB ×2, ccmC, ccmFC, ccmFN ×2). Three cytochrome c oxidase genes (cox1, cox2, cox3) and one maturase gene (matR) were also identified.

The variable gene set consists of three large ribosomal subunit genes (rpl2, rpl5, rpl16) and five small ribosomal subunit genes (rps3, rps4, rps7, rps12, rps14). Notably, rps12, rps14, and rpl16 are variably retained in angiosperm mitogenomes and are often considered non-canonical or plastid-derived, suggesting possible intracellular gene transfer or pseudogenization. The genome also encodes sdh4, a succinate dehydrogenase subunit, and mttB, a putative membrane transporter gene homologous to bacterial tatB, which may represent a case of horizontal gene transfer. The rRNA gene complement includes rrn5, rrn18, and rrn26, which constitute the core structural components of mitochondrial ribosomes. The tRNA set contains 18 distinct species, with multiple copies observed for trnE-UUC (×2), trnK-UUU (×2), trnM-CAU (×3), trnS-UGA (×2), and trnY-GUA (×4), providing redundancy typical of plant mitochondrial translation systems.

The mitochondrial PCGs of I indigotica display discernible codon usage biases, providing insights into the evolutionary dynamics and translational regulation within the mitogenome. A total of 61 sense codons encode the standard 20 amino acids, with the exception of methionine and tryptophan, both of which are encoded by single codons ( Figure 3A ). The distribution of relative synonymous codon usage (RSCU) values indicates varying degrees of codon preference, likely shaped by a combination of mutational biases, selective constraints, and tRNA availability.

Codons with RSCU values greater than 1 are considered preferentially used. Examples include GCT (alanine, 1.63), CGT (arginine, 1.26), AGA (arginine, 1.41), AAT (asparagine, 1.37), GAT (aspartic acid, 1.36), and TGT (cysteine, 1.24), suggesting favored usage during mitochondrial translation. In contrast, codons such as GCG (alanine, 0.49), CGG (arginine, 0.69), CGC (arginine, 0.53), AAC (asparagine, 0.63), and GAC (aspartic acid, 0.64) are underrepresented, with RSCU values below 0.7 ( Supplementary Table S2 ). These patterns collectively suggest a moderate codon usage bias in I. indigotica, with a tendency toward A- and T-ending codons at the third position, consistent with observations in other angiosperm mitogenomes.

A comprehensive analysis of SSRs in the mitogenome of I. indigotica revealed 51 distinct SSR loci, with repeat unit lengths ranging from 10 to 1,451 base pairs ( Figure 3B ). These SSRs were categorized into seven types—mononucleotide (p1), dinucleotide (p2), trinucleotide (p3), tetranucleotide (p4), pentanucleotide (p5), hexanucleotide (p6), and compound repeats (c)—and further grouped by repeat unit numbers (three to >7). Among them, compound SSRs were the most abundant (n = 14), followed by tetranucleotide repeats (n = 13) and dinucleotide repeats (n = 8) ( Supplementary Table S3 ). Mononucleotide and trinucleotide motifs were observed six times each, while pentanucleotide and hexanucleotide repeats were less frequent, with only two and one instance(s), respectively. The high prevalence of compound and tetranucleotide SSRs suggests their potential role in promoting genomic plasticity and structural variability. These SSRs were visualized by repeat type and length category in a 3D bar plot.

In addition to SSRs, other repetitive elements were also analyzed, including dispersed and tandem repeats ( Supplementary Table S4 ). A total of 396 dispersed repeats were identified, comprising 216 forward repeats and 180 palindromic repeats, with repeat lengths ranging from 30 bp to 11,396 bp. Tandem repeats were also detected, though less abundant. The spatial organization of these three major repeat types—SSRs (outermost circle), tandem repeats (middle circle), and dispersed repeats (innermost circle)—is illustrated in Figure 3C . In this circular map, green bars indicate SSR loci, light blue boxes represent tandem repeats, and arcs connect dispersed repeat pairs, with pink lines denoting forward repeats and purple lines indicating palindromic repeats. This layered distribution provides insight into the repeat architecture and potential recombinogenic regions across the I. indigotica mitogenome.

To experimentally validate homologous recombination events facilitated by repetitive elements, five repeat pairs—two forward repeats (F96 and F146) and three palindromic repeats (P93, P134, and P326)—were selected for PCR-based analysis ( Supplementary Table S5 ). Primers were designed to flank each repeat with ~100 bp of adjacent sequence. Two sets of primer combinations (F1/R1 and F2/R2) targeted the native configurations of each repeat region, while cross-combinations (F1/R2 and F2/R1) were designed to detect recombined conformations ( Table 2 ).

As shown in Figure 4 , agarose gel electrophoresis revealed distinct banding patterns for both native and recombinant configurations in multiple repeat pairs. For each repeat, at least two structural forms were detected: canonical conformations (F1R1 and F2R2) and recombined products (F1R2 and F2R1). The presence of clear bands in lanes corresponding to recombined primer pairs supports the hypothesis that repeat-mediated homologous recombination contributes to structural diversity in the I. indigotica mitogenome. Negative control lanes (NC) confirmed the specificity of the amplification reactions and excluded contamination. These findings highlight the dynamic nature of plant mitogenomes and the role of dispersed and palindromic repeats in generating structural variation.

The cpgenome of I. indigotica was successfully assembled with a total length of 153,821 bp and a typical quadripartite structure, comprising a large single-copy (LSC) region (83,577 bp), a small single-copy (SSC) region (17,702 bp), and two identical inverted repeat (IR) regions (26,271 bp each) ( Figure 2A ). Genome annotation identified gene clusters associated with photosystem I and II, ATP synthase, ribosomal proteins, and tRNAs, highlighting the conserved functionality of the cpgenome in photosynthesis and plastid gene expression. Functional categories were color-coded to distinguish among coding regions such as RNA polymerase subunits, NADH dehydrogenase components, and ribosomal RNA genes.

Comparative genome alignment revealed high conservation between the newly assembled cpgenome of I. indigotica and a previously reported sequence from the same species (OP620952), differing by only 6 base pairs in total length ( Figure 2B ). In contrast, alignment with the cpgenome of I. tinctoria (NC_028415) revealed minor sequence variations, suggesting evolutionary divergence between these congeneric species ( Figure 2C ). Further structural comparison of the LSC, SSC, and IR boundary regions between I. indigotica and I. tinctoria uncovered shifts in the positions of junctions such as JLB, JSB, JSA, and JLA, along with differences in adjacent gene lengths and arrangements ( Figure 2D ). These structural differences at IR boundaries, including variable placement of genes such as ycf1, rps19, and ndhF, underscore the dynamic nature of cpgenome evolution within the Isatis genus.

Intracellular DNA transfer is a common feature of plant genomes, contributing to genomic complexity and evolutionary adaptation. In the mitogenome of I. indigotica, a total of 17 mitochondrial plastid DNA transfers (MTPTs) were identified, with sequence lengths ranging from 74 bp to 1,367 bp and sequence identity between 74.4% and 99.8% ( Supplementary Table S6 ). These MTPTs collectively span several kilobases and reflect frequent plastid-to-mitochondrion DNA movement ( Figure 5A ). Among them, several MTPTs contained fragments of PCGs, such as rbcL (MTPT1), psaB (MTPT2, MTPT11), and ycf1 (MTPT3, MTPT4, MTPT9, MTPT10), indicating historical transfers of functional genetic elements. Notably, MTPT3 and MTPT4 represent bidirectional transfers of the same ycf1 fragment, with nearly complete sequence identity (99.8%) and reverse orientation.

In addition to PCGs, various ribosomal RNA and tRNA genes were also identified among MTPTs. For example, MTPT5 and MTPT6 encode partial sequences of rrn23, while MTPT7 and MTPT8 harbor homologs of rrn16, though with lower identity (~74%) and higher mismatch counts, suggesting older or less conserved transfers. tRNA-related MTPTs include trnP-UGG (MTPT12), trnW-CCA (MTPT13), trnD-GUC (MTPT14), and trnM-CAU (MTPT15, MTPT16), as well as trnI-CAU (MTPT17), showing varying levels of conservation. These findings indicate that not only protein-coding but also structural RNA elements were subject to intracellular translocation. While some tRNAs may retain functional integrity, others, particularly those with high mismatch and gap counts, may be relics of ancient transfers with limited biological function.

Beyond organelle–organelle transfer, comparative analysis based on the I. indigotica reference nuclear genome (NCBI accession VHIU00000000) revealed extensive integration of mitochondrial DNA (NUMTs) into the nuclear genome. As illustrated in Figure 5B , mitochondrial fragments were mapped to all seven nuclear chromosomes (chr1–chr7), with chr1, chr2, and chr5 exhibiting particularly dense insertion patterns. These NUMTs contribute to nuclear genome expansion and serve as evolutionary footprints of organellar-nuclear interaction. Together, the presence of MTPTs and NUMTs highlights the dynamic and ongoing exchange of genetic material between organellar and nuclear compartments in I. indigotica, providing important insights into organellar genome evolution and intracellular genomic plasticity.

The prediction-based strategy identified 380 putative C-to-U RNA-editing sites across mitochondrial PCGs of I.indigotica, based on conserved editing patterns from related species using PREPACT3 ( Supplementary Table S7 ). These predicted edits predominantly resulted in nonsynonymous changes, with a strong preference for specific amino acid substitutions. Notably, conversions from serine (S) to leucine (L), proline (P) to leucine (L), and serine (S) to phenylalanine (F) were the most frequent, suggesting potential functional shifts in hydrophobicity and structural properties of mitochondrial proteins. Gene-wise analysis revealed that nad4, ccmFn, and ccmB harbored the greatest number of editing sites, indicating that genes involved in respiratory and cytochrome c biogenesis processes are frequent targets of RNA editing ( Supplementary Figure S2 ). A detailed classification of RNA-editing events based on amino acid property transitions ( Supplementary Table S8 ) further highlights these tendencies. Of the 357 identified nonsynonymous edits, nearly half (48.74%) were transitions from hydrophilic to hydrophobic residues, indicating a substantial shift toward increased protein hydrophobicity. This group was dominated by conversions such as TCA (S) → TTA (L), TCT (S) → TTT (F), and CGG (R) → TGG (W). Hydrophobic-to-hydrophobic transitions accounted for 28.29%, mainly due to P → L and A → V substitutions. Hydrophilic-to-hydrophilic and hydrophobic-to-hydrophilic edits comprised 14.00% and 8.40%, respectively, with changes like CGT (R) → TGT (C) and CCT (P) → TCT (S) being prominent. Additionally, two editing events (0.56%) led to stop codons, potentially serving regulatory or degradation-related roles. This distribution emphasizes the functional importance of RNA editing in refining the hydrophobicity and polarity of mitochondrial proteins, thereby influencing their structural and biochemical dynamics.

In parallel, transcriptome read mapping against the I. indigotica mitogenome (C_AA108663) identified 230 RNA-editing sites ( Figures 5C, D ), derived from variant calling of publicly available RNA-seq data. Despite a lower total count, the editing spectrum overlapped substantially with prediction results. Key substitutions such as S→L, P→L, and S→F were consistently observed, and nad4 remained the most edited gene. These findings indicate that high-confidence editing events are reproducible across methodologies, and that a conserved subset of edits is likely to be biologically relevant. The concordant amino acid change patterns support the hypothesis that RNA editing modulates protein properties in a directionally selective manner, favoring alterations in polarity and hydrophobicity to maintain or optimize protein function within the mitochondrial environment.

Based on Figure 6 , the evolutionary and structural characteristics of the I. indigotica mitogenome were investigated through phylogenetic analysis, genome feature comparison, and multi-species alignment. The phylogenetic tree ( Figure 6A ), constructed from conserved mitochondrial PCGs across 26 representative angiosperm species, demonstrates that I. indigotica clusters tightly with I. tinctoria, confirming their close genetic relationship within the Brassicaceae family. Other genera, such as Raphanus, Brassica, and Lepidium, also form distinct, well-supported clades, reflecting the broader phylogenetic framework within the order Brassicales. The high bootstrap values across the tree nodes provide robust support for the inferred relationships and indicate strong conservation among mitochondrial genes across these taxa. This phylogenetic placement underscores the evolutionary consistency in mitogenomes within Brassicaceae and affirms the use of mitochondrial PCGs as reliable markers for resolving interspecies relationships. Furthermore, the comparison of mitogenome size and GC content among 18 Brassicaceae species ( Figure 6B ) reveals that although the total mitogenome size varies widely, ranging from approximately 220,000 to over 390,000 base pairs, the GC content remains relatively stable across species, typically clustering between 44% and 46%. Isatis tinctoria, marked in red, exhibits a genome size of approximately 283 kb with a GC content of ~45%, placing it near the average within its family. This trend suggests that while genome size can undergo lineage-specific expansion or reduction due to structural rearrangements or integration of foreign sequences, the GC content is more evolutionarily conserved and may be subject to functional constraints related to DNA stability or transcription efficiency.

In addition to phylogenetic and compositional features, the circular genome visualization ( Figure 6C ) offers a comprehensive view of the I. indigotica mitogenome and its structural architecture. The annotated internal rings display key genomic components, including coding sequences (CDS), transfer RNAs (tRNAs), ribosomal RNAs (rRNAs), GC content, and GC skew, which collectively reflect the functional landscape of the mitogenome. Surrounding these are twelve concentric outer rings that depict BLAST-based alignments between I. indigotica and twelve other selected species within Brassicaceae, enabling a comparative analysis of genome conservation and divergence. The presence of numerous conserved syntenic blocks indicates that many mitochondrial regions—particularly those associated with essential genes such as cox, nad, atp, and ccm—are strongly preserved across species, consistent with their critical roles in mitochondrial respiration and metabolism. At the same time, breaks and shifts in alignment patterns highlight lineage-specific rearrangements and structural variation, which may result from recombination, gene loss, or horizontal gene transfer events. These structural differences, although variable, often spare the core mitochondrial gene set, reinforcing the notion that while non-coding and intergenic regions are more evolutionarily plastic, the coding content remains stable to ensure essential mitochondrial functions. Taken together, this integrated analysis affirms that I. indigotica possesses a structurally stable mitogenome that shares extensive sequence similarity with its close relatives, while also harboring unique features that may reflect its individual evolutionary history. These findings provide valuable insight into mitogenome evolution within Brassicaceae and set a foundation for future investigations into mitochondrial gene function, inheritance, and adaptation.

The synteny and collinearity analyses of the mitogenome reveal notable patterns of both conservation and structural reconfiguration across different plant species. The syntenic comparisons, as illustrated in Figure 7A , show that extensive homologous sequence regions are retained between the focal genome and those of four other species: I. tinctoria, A. graveolens, T. mongolicum, and V. angularis ( Supplementary Figure S9 ). These homologous segments are represented by connecting curves, varying in intensity and continuity, which suggest differing degrees of sequence similarity and genome structural stability. The greatest density and continuity of these curves are observed between the focal genome and I. tinctoria, indicating a high degree of conservation in sequence composition and segmental arrangement. By contrast, the connections with the other three species appear more sporadic and fragmented, suggesting a comparatively lower level of sequence conservation and a higher degree of rearrangement. This discrepancy likely reflects differing evolutionary distances, with more distantly related species exhibiting greater divergence in mitogenome architecture.

Further insights into structural organization are provided by the collinearity analysis in Figure 7B , which focuses on three species: R. sativus, I. tinctoria, and the focal taxon. The visualization of LCBs, each demarcated by uniquely colored segments, facilitates a direct comparison of gene order, orientation, and overall genomic architecture. A pronounced similarity in block arrangement and orientation is evident between the focal genome and I. tinctoria, reinforcing the observation of strong genomic conservation inferred from the synteny analysis. In contrast, the alignment with R. sativus reveals a more disrupted block structure, characterized by numerous inversions, translocations, and possible insertions or deletions. These structural rearrangements result in a visually more complex and less orderly pattern of alignment, reflecting a greater divergence in genome structure. The differential collinearity among these three species suggests that mitogenome evolution in plants involves a balance between the retention of conserved functional elements and lineage-specific rearrangements. Such rearrangements may occur through recombination, duplication, or other mechanisms of genome plasticity, which contribute to the observed diversity in mitogenome organization across plant lineages.

To investigate the evolutionary dynamics of PCGs among angiosperms, three complementary analyses were conducted: gene copy number variation, Ka/Ks ratios, and Pi. Together, these results elucidate the extent of functional conservation, evolutionary constraint, and sequence variability across diverse lineages, with particular attention to the Brassicales and Fabales clades.

Gene copy number profiling ( Figure 8A ) revealed a heterogeneous landscape across the surveyed taxa. While the majority of PCGs were retained as single-copy genes in most species—particularly those involved in essential energy metabolism pathways such as ATP synthesis (e.g., atp1, atp6, atp9) and cytochrome complexes (e.g., cob, cox1–3)—a subset of genes exhibited lineage-specific duplications or absences. For instance, duplications in genes such as rps12, nad4L, and ccmFN were observed in several Brassicaceae members, suggesting possible subfunctionalization or compensatory mechanisms in mitochondrial translation and redox regulation. Meanwhile, frequent losses or underrepresentation of ribosomal protein genes (e.g., rpl2, rps1, rps14) in multiple lineages may indicate ongoing nuclear transfer or pseudogenization events. These patterns reflect both the functional rigidity and evolutionary plasticity of plant mitogenomes, where essential bioenergetic functions are maintained under strict constraint, while other genes tolerate structural fluidity.

The evolutionary constraints acting upon these genes were further evaluated through Ka/Ks analysis ( Figure 8B ), focusing on 16 representative Brassicaceae species. As expected for mitogenomes, purifying selection dominated across most PCGs, with Ka/Ks ratios predominantly below 1.0. However, exceptions were noted in genes such as ccmB, matR, and rps3, where elevated Ka/Ks values in specific taxa suggested relaxed selection or episodic positive selection. Notably, the variability observed in cytochrome c maturation genes (ccmB, ccmC, ccmFN) and ribosomal proteins (rps3, rps4) implies that these loci may be subject to adaptive pressures, potentially in response to environmental cues or cytonuclear co-evolution. Such divergence highlights regions of the mitogenome that may play a disproportionate role in lineage-specific adaptation or reproductive isolation.

Complementary to these findings, Pi analysis ( Figure 8C ) offered a broader population-level perspective on genetic variability. A clear disparity emerged between Brassicales and Fabales, with the latter exhibiting consistently higher Pi values across nearly all mitochondrial genes. Genes such as rps12, nad6, and ccmFC in Fabales showed particularly elevated polymorphism, indicating more rapid sequence turnover or relaxed constraints within this order. In contrast, Brassicales maintained low Pi, underscoring the genomic stability and strong conservation of mitogenomic content in this lineage. This contrast between lineages not only mirrors their differing evolutionary histories but also supports the notion that mitogenomes, though structurally conserved, can exhibit considerable heterogeneity in their molecular evolution.