In this study, young shoots with one bud and one leaf from the Camellia sinensis cv. ‘Zhuyeqi’, cultivated at the standardized tea planting demonstration base of Taohuayuan Agricultural Development Co., Ltd., Xinhua County, Hunan Province, China (27°47′N, 110°51′E), were used as experimental materials ( Figure 1 ). Genomic deoxyribonucleic acid (DNA) was extracted using the TIANGEN Plant Genomic DNA Extraction Kit (DP305), following the manufacturer’s standard protocol. DNA purity was assessed using a NanoDrop spectrophotometer (Thermo Scientific), with A260/A280 ratios between 1.8 and 2.0 and A260/A230 ratios ≥2.0. DNA integrity was verified using 1% agarose gel electrophoresis (5 V/cm, 30 min), ensuring DNA fragment lengths >30 kb to meet the high molecular weight DNA requirements for long-read sequencing (Xu et al., 2018; Zhang et al., 2019; Chen et al., 2023).

A combined strategy using second-generation (Illumina NovaSeq 6000) and third-generation (Oxford Nanopore PromethION) sequencing technologies was employed to balance read coverage and long-fragment assembly accuracy (Zhang et al., 2019; Wang et al., 2023). Second-generation sequencing library preparation: Genomic DNA was fragmented into 300–500 bp fragments using a Covaris system, followed by end repair, 3’ poly-A tailing, and Illumina adapter ligation. The library quality was assessed using a Qubit 3.0 fluorometer and an Agilent 2100 Bioanalyzer. Paired-end sequencing (2×150 bp) was performed on the NovaSeq 6000 platform (Zhang et al., 2019). Raw sequencing data were filtered using fastp (v0.23.4) to remove low-quality reads based on the following criteria:1) Reads containing adapter or primer sequences. 2) Reads with an average Phred quality score<Q5. 3) Reads with an N base proportion >5% (An et al., 2021; Li et al., 2024). Third-generation sequencing and quality control: DNA fragments >10 kb were enriched using a magnetic bead-based selection method, followed by end repair and adapter ligation using the SQK-LSK109 library preparation kit. Third-generation sequencing data were filtered using filtlong (v0.2.1) with parameters: –min_length 1000 and –min_mean_q 7. Perl scripts were used to analyze the read N50, average quality scores, and guanine-cytosine (GC) content to ensure the data met the requirements for high-accuracy assembly (Li et al., 2023b; Wang et al., 2023).

A hierarchical iterative assembly framework was constructed based on the high conservation of core plant mitochondrial genes (e.g., nad5, atp9, cox2) (Li et al., 2024). Third-generation sequencing data (PacBio/Nanopore) were initially aligned using minimap2 (v2.1). High-confidence seed sequences were selected based on alignment lengths >50 bp and coverage of at least three core genes to avoid interference from repetitive sequences (Li et al., 2024). Seed sequences were extended through multiple rounds of iterative alignment (overlap >1 kb) and integrated into an initial assembly dataset. To address the unique subcircular and non-circular topologies of tea tree mitochondria (Zhang et al., 2019; Bi et al., 2022; Li et al., 2024), third-generation sequencing data were corrected using Canu, and a hybrid assembly strategy was implemented by combining second-generation data (mapped using Bowtie2 v2.3.5.1) with Unicycler (v0.4.8). This approach significantly improved haplotype coverage depth and assembly continuity (Wang et al., 2024b). Further analysis using the visualization tool Bandage (v0.8.1) was performed to examine contig branching patterns. Manual corrections were applied to refine multichromosomal structure boundaries, based on differences in coverage depth and the distribution of repeat sequences. This process enabled the precise reconstruction of the complex multichromosomal configuration (Bi et al., 2022; Li et al., 2024). This strategy effectively overcame the fragmentation issues caused by high repeat sequences, providing a solid structural foundation for subsequent adaptive evolution analyses (Li et al., 2024; Wang et al., 2024b).

A multidimensional annotation strategy was applied: (1) Core gene annotation: The mitochondrial genomes of closely related species were used as a reference database. Homologous sequences were identified using BLAST (E-value<1e-5). Exon-intron boundaries were optimized using Geneious Prime, and chloroplast homologous fragments were manually excluded (threshold: sequence similarity ≥80% (Lu et al., 2022; Li et al., 2023b). (2) Transfer RNA (tRNA) identification: tRNAs were predicted using tRNAscan-SE (v2.0) under the default eukaryotic mode, with secondary validation by ARWEN to ensure stability of tRNA secondary structures (Zhang et al., 2019; Li et al., 2024). (3) Open reading frame (ORF) screening and functional annotation: Non-redundant sequences ≥102 bp were identified using ORF Finder (National Center for Biotechnology Information, NCBI). Overlapping regions with known genes and low-complexity regions (analyzed with RepeatMasker) were excluded. ORFs ≥300 bp were annotated using BLASTX against the nr database, focusing on gene clusters related to adaptive evolution Complete chloroplast genome sequence of Camellia sinensis genome structure, adaptive evolution, and phylogenetic relationships (Li et al., 2024).

To systematically analyze the functional characteristics of mitochondrial RNA editing sites in the Camellia sinensis cv. ‘Zhuyeqi,’ this study employed an integrated multi-omics strategy. First, RNA-seq data generated from the Illumina platform were aligned to mitochondrial genome coding regions using Bowtie2 (v2.3.5.1) with high precision. The alignment results in BAM format were then sorted and filtered for quality using SAMtools (v1.9) (Li et al., 2009; Danecek et al., 2021). Next, single nucleotide polymorphism (SNP) sites between genomic DNA and transcriptome data were identified using BCFtools (v1.9-170). These results were integrated with predictions from the Predictive RNA Editor for Plant Mitochondria to construct a high-confidence set of RNA editing sites (Mower, 2009; Danecek et al., 2021). To further validate the biological significance of editing events, the Kyte-Doolittle hydropathy index was used to quantify changes in the physicochemical properties of amino acids induced by editing (Zhang et al., 2022). Finally, a generalized linear model was applied to analyze the association between RNA editing sites and signals of positive selection (Ka/Ks >1), revealing the regulatory role of RNA editing in the adaptive evolution of tea tree mitochondria (Zhang et al., 2022; Li et al., 2023b).

To investigate the codon usage bias and the evolutionary driving forces in the mitochondrial genome of the tea plant variety ‘Zhuyeqi’, this study used the Relative Synonymous Codon Usage (RSCU) index for quantitative analysis. The RSCU is calculated using Equation 1:

Here, f o b s represents the observed frequency of a synonymous codon, and f e x p     represents the expected frequency under a neutral evolutionary model. An RSCU value of 1 indicates neutral selection, while RSCU > 1 or< 1 suggests positive or negative selection pressure, respectively (Zhang et al., 2022; Li et al., 2023b). To eliminate redundancy from homologous sequences, unique mitochondrial coding sequences were filtered using CD-HIT (identity = 100%). The RSCU values were then calculated and statistically tested in batches using a custom Perl script, available at https://github.com/ari-dasci/S-TSFE-DL. The expected frequencies were based on a uniform distribution model (Fu et al., 2012).

To systematically analyze the characteristics of repetitive sequences in the mitochondrial genome, this study utilized a multi-dimensional strategy for repetitive sequence detection: (1) Simple Sequence Repeats (SSR): Whole-genome scanning for SSRs was performed using the MISA software (v1.0). The parameters were set based on plant mitochondrial genome research standards. The minimum repeat numbers for mono- to hexanucleotide repeat units were set as 10, 5, 4, 3, 3, and 3, respectively, to balance detection sensitivity and specificity (Sahu et al., 2012; Zhang et al., 2019). (2) Tandem Repeats: Tandem repeats were analyzed using the TRF algorithm (v4.0.9). The parameters were set as follows: match score = 2, mismatch penalty = 7, indel penalty = 7, and minimum alignment score = 80. This ensured the effective identification of periodic repeat units ≥50 bp in length (Benson, 1999; Wynn and Christensen, 2019). (3) Dispersed Repeats: Dispersed repeats were detected by performing a whole-genome self-comparison using BLASTn (v2.10.1), with optimized settings for word size (word_size = 7) and significance threshold (evalue = 1e-5). The REPuter algorithm was then used to remove redundant signals and to exclude tandem repeat interference (Wynn and Christensen, 2019; Liu et al., 2024). Finally, a multi-level interaction network map of repetitive sequences was constructed using Circos (v0.69-5), revealing the spatial distribution patterns of repetitive elements within the genome (Zhang et al., 2019).

Based on the homologous gene clusters of the following reference species: Zea mays subsp. mays (NC_007982.1), Brassica oleracea (NC_016118.1), Gossypium barbadense (NC_028254.1), Gossypium arboreum (NC_035073.1), Styphnolobium japonicum (NC_039596.1), Chrysanthemum boreale (NC_039757.1), and Camellia sinensis (NC_043914.1), an evolutionary pressure model was constructed. First, global multiple sequence alignment was performed using MAFFT v7.427 (Katoh and Standley, 2013), and low-conserved regions were removed with Gblocks. Then, Ka and Ks values were calculated using the MLWL algorithm in KaKs_Calculator v2.0 (Li et al., 2024), with the significance threshold set at p< 0.05. Particular focus was given to genes under positive selection (Ka/Ks > 1) and genes under strong purifying selection (Ka/Ks< 0.5), which were further analyzed for functional associations with previously reported core functional genes in tea plant mitochondria (Li et al., 2024). Boxplots of Ka/Ks values across evolutionary branches were constructed using the ggpubr package in R, and an evolutionary model of selection pressure over time and across lineages was integrated with PhyloTree (Wang et al., 2024b).

Based on the MAFFT multiple sequence alignment (parameters: –retree 1 –maxiterate 1000) of homologous gene clusters from the following reference species: Zea mays subsp. mays (NC_007982.1), Brassica oleracea (NC_016118.1), Gossypium barbadense (NC_028254.1), Gossypium arboreum (NC_035073.1), Styphnolobium japonicum (NC_039596.1), Chrysanthemum boreale (NC_039757.1), and Camellia sinensis (NC_043914.1), nucleotide polymorphism analysis was performed using the sliding window method (window length: 500 bp, step size: 100 bp) (Li et al., 2024). The average Pi value of single nucleotide polymorphism sites was calculated using DnaSP v5.0, and confidence intervals were assessed through 1,000 bootstrap resampling iterations to quantify the level of genetic differentiation among species (Li et al., 2024; Wang et al., 2024b). To eliminate interference from low-quality alignment regions, highly variable sites were filtered prior to alignment using Gblocks v0.91b (allowing gaps with the parameter: -b5=h) (Wang et al., 2024b).

Based on the multiple sequence alignment of whole-genome coding sequences generated by MAFFT, conserved regions were selected using trimAl v1.4 (parameters: -gt 0.7 -cons 60) (Li et al., 2024) to retain high-confidence sites. The optimal evolutionary model was determined using PartitionFinder v2.1.1, supporting the GTR+GAMMA+I model (Lanfear et al., 2016). A maximum likelihood tree was constructed with RAxML v8.2.10 (parameters: -m GTRGAMMAI -f a -x 12345), and Bayesian inference was performed using MrBayes v3.2.6 (number of chains = 4, sampling frequency = 1,000, convergence criterion: average standard deviation< 0.01). Finally, 1,000 bootstrap support values and posterior probabilities were integrated to enhance the reliability of the tree topology (Li et al., 2024; Wang et al., 2024b).

Whole-genome synteny alignment was performed using nucmer v4.0 (parameters: –maxmatch -c 65 -l 40), and syntenic blocks and rearrangement events were identified with SyRI v1.6 (Li et al., 2023b). High-resolution dot plots were generated using ggplot2 (Wickham, 2017). For mitochondrial-to-chloroplast horizontal transfer sequences, BLASTN (threshold: E-value ≤ 1e-10, coverage ≥ 80%) (Altschul et al., 1997) was combined with GeneWise v2.4.1 to verify open reading frame integrity (Zhang et al., 2019) and exclude false-positive transfer events. Cross-genome interaction networks were visualized with CIRCOS v0.69-5, and the dynamics of genome structure were revealed based on SSRs and long terminal repeats identified by REPuter v2.74 (Kurtz et al., 2001; Li et al., 2023b; Liu et al., 2024).

This study is the first to elucidate the complex multichromosomal configuration of the mitochondrial genome in the Camellia sinensis cv. ‘Zhuyeqi’ (GenBank accession numbers: PQ763484–PQ763490). The genome exhibits a unique structure composed of one circular chromosome and six linear chromosomes ( Figure 2 ; Supplementary Table 1 ), with a total length of 911,255 bp (individual chromosome lengths ranging from 133 to 311,104 bp) and a GC content of 46% ( Supplementary Table 2 ). The genome size and GC content are highly conserved compared to previously reported Camellia species (701,719–1,081,966 bp, GC 44–46%), supporting the evolutionary stability of mitochondrial genomes in the Theaceae family (Li et al., 2023b). Gene annotation identified 77 functional genes, including 38 protein-coding genes (PCGs), 33 tRNAs, 3 ribosomal RNAs (rRNAs), and 3 pseudogenes ( Table 1 ; Supplementary Table 3 ). The core functional modules include adenosine triphosphate (ATP) synthase (atp1/4/6/8/9), nicotinamide adenine dinucleotide (NADH) dehydrogenase (nad1–7/9), and cytochrome c oxidase (cox1–3), consistent with the conserved characteristics of mitochondrial genomes in terrestrial plants (Lu et al., 2022; Chen et al., 2023). Notably, the GC content of PCGs (43%) is significantly lower than that of tRNAs (50%) and rRNAs (48%) ( Supplementary Table 2 ). This nucleotide composition heterogeneity may be associated with functional differentiation or post-transcriptional regulatory mechanisms (Chen et al., 2023; Li et al., 2024).

Structural analysis of the genome revealed multiple copy numbers for some genes, such as trnM-CAT (5 copies), trnP-TGG (3 copies), and trnF-GAA, trnI-GAT, and four other genes with two copies each. These duplications suggest that recombination events mediated by repeat sequences or horizontal gene transfer might drive genome structural variation (Li et al., 2024; Liu et al., 2024). Additionally, the presence of the sdh3 pseudogene may reflect functional redundancy and differentiation of the mitochondrial genome during evolution, aligning with the adaptive evolutionary strategies of the Camellia genus (Li et al., 2023b, 2024). Moreover, the distribution of introns exhibits high heterogeneity: nine genes, including ccmFc and rpl2, contain a single intron; nad4 contains two introns, while nad1/2/5/7 each contain four introns. This complex splicing pattern may contribute to the fine regulation of gene expression through alternative splicing, closely related to the dynamic plasticity of the mitochondrial transcriptome in tea plants (Zhang et al., 2019, 2022). Similar intron structures have also been reported in Camellia sinensis var. assamica, suggesting that these may represent conserved regulatory features of mitochondrial genomes in the Theaceae family (Zhang et al., 2019).

This study systematically analyzed 38 PCGs in the mitochondrial genome of the Camellia sinensis cv. ‘Zhuyeqi’ and identified a total of 556 RNA editing sites ( Figure 3 ), with C-to-U editing being the predominant type. This result is consistent with the dominant C-to-U editing pattern observed in mitochondrial genomes of angiosperms (Li et al., 2021; Zhang et al., 2022). The editing sites were unevenly distributed among the genes: the ccmFn gene contained the most editing sites (38), followed by ccmB (34), mttB (33), and nad2 (30), while sdh3 had only 2 sites. Notably, the high density of editing sites in the ccmB gene (34 sites) aligns with reports from other tea cultivars (Zhang et al., 2022), and its features under positive selection may reflect its critical role in mitochondrial energy metabolism and its adaptive evolutionary demands (Li et al., 2024; Zhang et al., 2024). Additionally, the high number of editing sites in ccmFn is similar to the editing hotspots observed in mitochondrial genes of plants like Populus, suggesting that strong selective pressure due to functional conservation may act on this gene during cross-species evolution (Zhang et al., 2022, 2024).

The amino acid substitutions caused by RNA editing were classified into five categories: hydrophilic-to-hydrophilic, hydrophilic-to-hydrophobic, hydrophilic-to-stop codon, hydrophobic-to-hydrophilic, and hydrophobic-to-hydrophobic ( Table 2 ). Among these, hydrophilic-to-hydrophobic substitutions were the most frequent (48%), with typical examples including TCA (S, serine) → TTA (L, leucine), TCG (S) → TTG (L), and CGG (R, arginine) → TGG (W, tryptophan). These substitutions significantly enhanced the local hydrophobicity of proteins, which may contribute to remodeling transmembrane domains or the hydrophobic cores of enzyme active sites, thereby regulating the assembly efficiency and functionality of mitochondrial oxidative phosphorylation complexes (Li et al., 2021; Zhang et al., 2022). Hydrophobic-to-hydrophobic substitutions accounted for 31% (e.g., CCA → CTA causing P → L, GCC → GTC causing A → V). Although these substitutions did not change amino acid polarity, differences in side-chain volume and hydrophobicity may lead to subtle conformational adjustments, such as altering the packing density of ATP synthase transmembrane helices or the spatial constraints of substrate-binding pockets (Zhang et al., 2019; Duan and Lai, 2020). Hydrophilic-to-stop codon substitutions were the least frequent (1%), such as CAG (Q) → TAG (X) and CGA (R) → TGA (X).

This study conducted a codon usage analysis on 10,828 codons from 38 PCGs in the mitochondrial genome of the Camellia sinensis cv. ‘Zhuyeqi’. The results revealed a significant codon usage bias pattern. Among aliphatic amino acids, leucine (Leu) had the highest expression frequency (1,113 occurrences, accounting for 10%), followed by serine (Ser, 9%) and arginine (Arg, 7%), whereas tryptophan (Trp) showed a significantly lower usage frequency (1%) compared to other amino acids ( Figure 4 ; Supplementary Table 4 ). This distribution pattern is highly consistent with the conserved evolutionary trends of mitochondrial genomes in Camellia species (Mebeaselassie and Zhu, 2022; Li et al., 2024). Notably, among aromatic amino acids, the phenylalanine (Phe) codon UUU exhibited extreme preference (379 occurrences, 4%), with a 23% higher usage frequency than its synonymous codon UUC (3%). This observation aligns with the preferential selection of A/T-ending codons in the chloroplast genomes of Theaceae plants (Wang et al., 2022). Similarly, the isoleucine (Ile) codon AUU (3%) and the glutamic acid (Glu) codon GAA (3%) also displayed high expression frequencies. In contrast, no expression activity was detected for the methionine (Met) codons CUG/UUG, suggesting a potential functional deficit in the mitochondrial translation system (Mebeaselassie and Zhu, 2022; Chen et al., 2023).

Across the genome, 31 codons with relative synonymous codon usage (RSCU) values >1 were identified, accounting for 47% of all codons. Among these, the methionine start codon AUG had the highest RSCU value (3.0) in the genome. This result is consistent with the conserved functional characteristics of start codons reported by Li et al. (2024) in a comparative mitochondrial genomics study of tea variants (Li et al., 2024). Analysis of effective number of codons (ENC) and GC3s data ( Supplementary Table 5 ) showed that high-RSCU codons were predominantly found in AU-rich regions with low GC content (e.g., AUU, UUU). This is highly related to the adaptive evolutionary strategy of terrestrial plant mitochondrial genomes, which reduce translation energy costs through codon deoptimization (Wang et al., 2022; Li et al., 2024).

The mitochondrial genome of the Camellia sinensis cv. ‘Zhuyeqi’ exhibits a typical multilayered repetitive sequence structure, including SSRs, tandem repeats, and dispersed repeats. Quantitative analysis identified a total of 1,078 repetitive sequences, with SSRs (269, 25%) and tandem repeats (29, 3%) as minor components, while dispersed repeats (780, 72%) were the predominant type ( Figure 5 ; Supplementary Table 6 ). This distribution pattern aligns with the conserved characteristics of tea mitochondrial genomes, suggesting that dispersed repeats play a central role in driving the dynamic evolution of the mitochondrial genome (Li et al., 2024).

The total length of dispersed repeats was 68,976 bp, accounting for 8% of the genome, which is significantly higher than reported for other tea cultivars (e.g., C. sinensis var. assamica) (Huang et al., 2014; Li et al., 2023b). The distribution showed marked heterogeneity among chromosomes, with chr4 (8,082 bp), chr6 (3,900 bp), and chr1 (1,281/1,277 bp) having the largest alignment lengths. This may be associated with localized recombination suppression or active transposon activity (Huang et al., 2014). Subtype analysis revealed that forward repeats (354, 45%) and palindromic repeats (426, 55%) were the major types ( Supplementary Table 7 ). Notably, chr4 contained the largest forward repeat (8,082 bp), while chr1 carried the largest palindromic repeat (1,281 bp). These findings support the hypothesis that mitochondrial genomes maintain structural stability through homologous recombination (Zhang et al., 2019; Li et al., 2023b). Short dispersed repeats (<100 bp) accounted for 86% (668) of all dispersed repeats. Their high abundance suggests that DNA slippage replication may drive genome microevolution, a phenomenon widely observed in the adaptive evolution of plant mitochondrial genomes (Li et al., 2023a, 2024).

Satellite DNA analysis identified 29 tandem repeat units with lengths ranging from 6–39 bp, showing significant asymmetry in chromosome distribution: chr1 was enriched with 11 sites, chr3 had 5, chr6 only 1, and chr7–11 lacked tandem repeats entirely ( Supplementary Table 8 ). This uneven distribution is associated with the compartmentalized multichromosomal structure of mitochondria and may reflect the regulation of mitochondrial genome 3D conformation by nuclear-mitochondrial interactions (Huang et al., 2014; Liu et al., 2024).

In molecular evolution studies, the ratio of Ka to Ks is a key indicator for assessing the selective pressure on protein-coding genes. This study systematically compared the evolutionary dynamics of 37 core functional genes in the mitochondrial genome of the Camellia sinensis cv.’Zhuyeqi’with those of representative species, including Camellia sinensis (NC_043914.1), Brassica oleracea (NC_016118.1), Chrysanthemum boreale (NC_039757.1), Gossypium arboreum (NC_035073.1), G. barbadense (NC_028254.1), Styphnolobium japonicum (NC_039596.1), and Zea mays subsp. mays (NC_007982.1) ( Figure 6 ; Supplementary Table 9 ). The analysis revealed cross-lineage adaptive evolutionary characteristics.

The Ka/Ks values for the 37 PCGs ranged from 0.09 to 2.70 ( Supplementary Table 10 ), consistent with the Ka/Ks distribution patterns observed in the mitochondrial genome of Camellia species [2]. Notably, the atp9 gene of Zea mays subsp. mays showed the lowest Ka/Ks value (0.09), indicating strong purifying selection, likely associated with the functional conservation of the ATP synthase complex in energy metabolism (Xia et al., 2014; Huang et al., 2022).

In contrast, the sdh4 gene of Chrysanthemum boreale exhibited the highest Ka/Ks value (2.70), followed by the ccmB gene (2.22), suggesting that these genes may have undergone positive selection. This finding aligns with the positive selection observed for Nad1 and Sdh3 genes in the mitochondrial genome of Camellia species, highlighting the convergent evolution of key metabolic genes (e.g., succinate dehydrogenase complex and cytochrome c maturation-related genes) in plant adaptive evolution (Li et al., 2024). Notably, the average Ka/Ks value for the ccmB gene was 1.38, significantly exceeding the neutral evolution threshold (Ka/Ks = 1). Its positive selection characteristics appear conserved across multiple plant lineages (Li et al., 2023b, 2024).

This study conducted a phylogenetic analysis of nucleotide diversity (Pi) for 40 core functional genes in the mitochondrial genome of the Camellia sinensis cv. ‘Zhuyeqi’ ( Figure 7 ; Supplementary Table 11 ). The results revealed multidimensional heterogeneity in genetic variation. The Pi values of functional genes showed a wide range (0–0.09), with the ribosomal protein genes rps4 (Pi = 0.09), ATP synthase subunit gene atp8 (0.09), and rps1 (0.08) exhibiting exceptionally high levels of variation. In contrast, rps19 (Pi = 0) was completely conserved, suggesting it is under strong functional constraint or purifying selection (Ka/Ks< 1) (Li et al., 2024). Notably, the succinate dehydrogenase gene sdh3 (Pi > 0.07) and rps14 (Pi > 0.07) displayed high levels of variation, which is consistent with previous studies reporting positive selection (Ka/Ks > 1) in the nad1 and sdh3 genes (Li et al., 2024), further supporting the critical role of mitochondrial genomes in the adaptive evolution of tea plants (Li et al., 2024).

Among rRNA genes, rrn5 (Pi = 0.03) showed significantly higher variation compared to rrn26 (Pi = 0.02), while the absence of Pi values for rrn18 may be attributed to its highly conserved secondary structure and RNA editing mechanisms. Studies have shown that mitochondrial RNA editing is predominantly C-to-U, with a strong preference for T bases in the surrounding sequences (Li et al., 2024). This editing pattern likely contributes to maintaining functional stability by regulating translation efficiency (Zhang et al., 2019). Comparative analysis of genomic structures revealed significant differences in variation patterns between mitochondrial and chloroplast genomes. The chloroplast genome exhibited stronger conservation in inverted repeats (IRs) (Pi< 0.005) (Peng et al., 2021), whereas the mitochondrial genome showed increased structural plasticity due to frequent recombination and horizontal gene transfer events (Li et al., 2024).

This study systematically analyzed the homologous sequences between the mitochondrial and chloroplast genomes of the Camellia sinensis cv. ‘Zhuyeqi’. The mitochondrial genome, with a total length of 911,255 bp, is approximately 5.8 times larger than the chloroplast genome (157,072 bp) ( Figure 8 ; Supplementary Table 12 ), consistent with the reported size variation of tea mitochondrial genomes (701,719–1,081,966 bp) (Li et al., 2023b). Whole-genome alignment identified 66 homologous fragments between the chloroplast and mitochondrial genomes, ranging from 32 to 7,721 bp, with a total length of 25,656 bp ( Supplementary Table 13 ). These homologous sequences account for 15% of the chloroplast genome and 2% of the mitochondrial genome. The proportion of mitochondrial homologous sequences (2%) falls within the typical range (1.91%-2.45%) reported for other tea cultivars (Li et al., 2024). This proportion differs significantly from legumes, such as Dalbergia odorifera, where chloroplast-derived sequences account for approximately 4% of the mitochondrial genome (Hong et al., 2021), but both systems highlight the widespread occurrence of inter-organellar DNA transfer (Wang et al., 2012; Li et al., 2023b). It is important to note that the evolutionary rates and recombination mechanisms of plant mitochondrial and chloroplast genomes differ significantly (Sloan et al., 2012; Tyszka et al., 2023), and no direct quantitative relationship has been established between the proportion of homologous sequences and the activity of DNA transfer (Li et al., 2024; Wang et al., 2024a). Therefore, the observed proportion of homologous sequences may represent a baseline level of gene flow between organelles in tea plants, while the activity of the specific transfer mechanisms requires further analysis, including recombination frequency and insertion site characteristics (Hong et al., 2021; Li et al., 2024).

In the chloroplast genome, homologous sequences (23,747 bp, 15%) contain 27 complete genes, including key functional genes such as trnA-UGC, rrn16, and rpl2. Additionally, the homologous regions of rps12, atpE, and rrn16 covered 75%, 91%, and 58% of their original gene lengths, respectively ( Supplementary Table 13 ). This suggests that chloroplast DNA transferred to the mitochondria may preferentially retain functional domains (Li et al., 2023b). In contrast, the mitochondrial homologous sequences (14,776 bp, 2%) included only 14 complete genes (e.g., trnV-GAC and trnW-CCA), with some homologous genes showing coverage below 20% ( Supplementary Table 13 ). This indicates that the mitochondrial genome may impose stricter selection or degradation mechanisms on transferred sequences (Li et al., 2024). Such differences may be driven by the strong purifying selection acting on the mitochondrial genome (Li et al., 2024).

A maximum likelihood phylogenetic tree constructed based on 38 mitochondrial protein-coding genes ( Figure 9 ) revealed five monophyletic groups among 26 representative species. The study selected three Poaceae species, Triticum timopheevii (NC_022714.1), Chrysopogon zizanioides (NC_0656367.1), and Zea mays (NC_007982.1), as outgroups (Yang et al., 2024a). The remaining species clustered into branches corresponding to Theaceae, Asteraceae, Brassicaceae, Malvaceae, and Fabaceae. Notably, the Camellia sinensis cv. ‘Zhuyeqi’ clustered wit1h 100% bootstrap support within the Camellia branch, forming a monophyletic group with C. sinensis (OM809792.1) and the Assam variety C. sinensis var. assamica (OL989850.1). This result aligns with the monophyly of the Camellia tribe observed in phylogenetic analyses based on the matR gene (Yang et al., 2006). Genetic differentiation distance analysis showed that Camellia sinensis cv. ‘Zhuyeqi’ shares a closer relationship with three major cultivated tea varieties—C. sinensis var. pubilimba (ON782577.1), C. sinensis (NC_043914.1), and C. sinensis (MH376284.1)—than the inter-family differentiation threshold. This finding suggests that these varieties may have originated from a shared ancestral population that underwent rapid radiation evolution (Xu, 2019; Li et al., 2024).

This study conducted a synteny analysis of the mitochondrial genome of the Camellia sinensis cv. ‘Zhuyeqi’ using whole-genome alignment strategies ( Figure 10 ; Supplementary Table 14 ). The results revealed significant structural conservation between Camellia sinensis cv. ‘Zhuyeqi’ and other Camellia sinensis species, with homologous regions covering 64% (579,821 bp) of the Camellia sinensis cv. ‘Zhuyeqi’ genome and 88% (620,962 bp) of the other Camellia genome. This observation contrasts sharply with the frequent gene rearrangements and heterogeneous repeat sequences commonly found in Camellia mitochondrial genomes (Lu et al., 2022; Li et al., 2024). Notably, the structural conservation in the Camellia sinensis cv. ‘Zhuyeqi’ may be attributed to artificial selection pressures during domestication, which preserved core functional modules. This aligns closely with the “domestication bottleneck effect” theory proposed by Li et al. (2024) (Li et al., 2024). Specifically, although the mitochondrial genome of Camellia species can contain up to 44% SSRs and exhibit numerous re6combination hotspots (Li et al., 2023b), the observed conservation in this study may result from the selective stabilization of key functional units such as respiratory chain complexes (Wei et al., 2018). Cross-species comparisons showed that syntenic regions between Camellia sinensis cv. ‘Zhuyeqi’ and the monocot Zea mays subsp. mays accounted for only 9% (77,430 bp) of the genome. This difference is likely due to genome remodeling events triggered by early divergences among angiosperms (Wei et al., 2018; Yang et al., 2024a).