The A. cerana-Southern China samples were obtained from the Hongchang Town, Chaonan District, Shantou, Guangdong Province, China (N 23.31, E 115.95). Fifty drone pupae were collected from one colony, then placed in absolute ethanol and stored in a −20 °C freezer. The genomic DNA was extracted and sent to Berry Genomics Co., Ltd. (Beijing, China) for PacBio HiFi sequencing.

PacBio HiFi reads generated for A. cerana–Southern China were first used to assemble the nuclear genome, including mitochondrial sequences, using Hifiasm v0.19.5-r587 (Cheng et al., 2021; Cheng et al., 2022). The resulting contigs from the Hifiasm assembly were indexed to construct a local BLAST database. A reference A. cerana mitochondrial genome (GenBank accession NC_014295) was downloaded and used solely to identify mitochondrial contigs within the Hifiasm assembly by BLAST searches, rather than for genome assembly. The identified mitochondrial contigs were then extracted for downstream analyses. Preliminary annotation of the mitochondrial genome was performed using the MITOS Web Server (Bernt et al., 2013). The annotation results were subsequently refined and validated using GeSeq implemented in Chlorobox (Tillich et al., 2017) by comparison with previously reported protein-coding and RNA genes from related species. A circular map of the mitochondrial genome was generated using Chloroplot (Zheng et al., 2020). The secondary structures of tRNA genes were predicted using the MITOS Web Server under the invertebrate mitochondrial genetic code (Bernt et al., 2013), and the tRNA secondary structure diagrams were manually redrawn using Adobe Illustrator.

BioEdit (version 7.0.9.0) (Hall, 1999) was used to calculate the base composition and codon usage of the complete A. cerana-Southern China mitochondrial genome sequence and the protein-coding genes (PCGs), tRNA genes, rRNA genes and AT-rich regions. Codon usage bias, represented by relative synonymous codon usage (RSCU), was calculated using codonW software (retrieved from https://codonw.sourceforge.net/) and was plotted into a bar graph using JSHYCloud (http://cloud.genepioneer.com:9929). The RSCU greater than 1 indicated a preference for certain amino acids. The AT and GC biases were calculated according to the following formulas: AT skew = (A-T)/(A + T) and GC skew = (G-C)/(G + C).

The mitochondrial genome of A. cerana–Southern China, assembled from PacBio HiFi sequencing data, has been deposited in the NCBI GenBank database under the accession number PP692293. In addition to this newly generated mitogenome, a total of 18 honey bee mitochondrial genomes were included in this study, comprising 15 previously published sequences retrieved from GenBank. The mitochondrial genomes of Apis mellifera sinisxinyuan and Apis mellifera capensis were selected as outgroups (Table 1). Because mitochondrial tRNA and rRNA genes are highly conserved in both nucleotide sequences and secondary structures, their inclusion provides limited phylogenetic resolution. Therefore, 13 mitochondrial protein-coding genes (PCGs) were used to infer the phylogenetic relationships among A. cerana populations from different geographic regions (Figure 1).

The 13 PCG sequences were input into PhyloSuite (V1.2.3) software (Xiang et al., 2023; Zhang et al., 2020) and aligned using MAFFT. Then, the PCG sequences were concatenated. Modelfind was used for evolutionary model selection with a BIC value. The evolutionary tree was constructed using MrBayes and an IQ tree with default parameters.

The PacBio HiFi–assembled nuclear genome was anchored to chromosomes using RagTag (Alonge et al., 2022) against the Apis cerana chromosome-level reference (GCA_029169275.1), after which the complete mitochondrial genome served as the BLASTn query to assess mito-nuclear collinearity (synteny).

Using PacBio HiFi sequencing data, we assembled a high-quality mitochondrial genome of A. cerana–Southern China (GenBank accession PP692293). PacBio HiFi sequencing generated 64.44 Gb of high-quality long-read data from drone pupae, with an average read length of approximately 17.8 kb and a mapping ratio of 99.90%, providing sufficient coverage for reliable genome assembly. The Hifiasm assembly produced a high-quality nuclear genome (total length: 231.9 Mb) with 40 contigs and a contig N50 of 11.78 Mb, supported by a BUSCO completeness score of 97.90%.

We assembled the complete A. cerana–Southern China (GenBank: PP692293) mitochondrial genome (Figure 2). The complete mitochondrial genome is 16,137 bp in length and contains 37 genes, including 13 protein-coding genes (PCGs), 22 tRNA genes, two rRNA genes, and an AT-rich region (Table 2). Among these genes, 23 are located on the J strand (majority strand), whereas the remaining 14 are encoded on the N strand (minority strand) (Figure 2; Table 2). The nucleotide composition of the complete mitochondrial genome is A = 42.41%, T = 41.82%, C = 9.59%, and G = 6.18%, resulting in a high A + T content of 84.23%. Most protein-coding genes (PCGs) initiated with standard ATN start codons, while several genes used TTA as an alternative initiation codon. In several genes, including ND5, ND4, ND4L, and ND1, TTA is annotated as an alternative initiation codon, based on conserved open reading frames and comparative alignment with closely related Apis mitochondrial genomes.

Regarding termination codons, most PCGs terminate with the conventional mitochondrial stop codon TAA. In a small number of cases, truncated stop codons were inferred, which are presumed to be completed through post-transcriptional polyadenylation, a common mechanism in insect mitochondrial genomes. The intergenic regions of the A. cerana–Southern China mitochondrial genome range from 0 to 229 bp in length.

The total length of the 13 PCGs in A. cerana-Southern China was 11,055 bp, accounting for 68.51% of the entire mitochondrial genome (Table 3). The length and base composition of the PCGs were varied, the 13 PCGs ranged from 162 bp (ATP8) to 1668 bp (ND5). Except for the A + T content of COX1 and COX2, the A + T content of the other PCGs was greater than 80%, indicating an AT preference. According to the AT-skew and GC-skew results, most of the PCGs were biased toward T, while all PCGs showed a preference for codon C. We conducted an analysis of the frequency and RSCU values for 13 PCGs in the A. cerana-Southern China mitochondrial genome. The RSCU values for the 13 PCGs were greater than 1, except for those of the two amino acids Met and Trp, which were each encoded by only one codon (Figure 2). This suggests that all codons were biased. The TTA codon had the highest RSCU value in A. cerana-Southern China (3.84), followed by the AGA codon (3.46). NNA-type codons accounted for 41.93% of the total number of codons, indicating that codons with an A at the third position were the most frequent. This study conducted a synteny analysis of the 13 protein-coding genes (PCGs) across 18 mitochondrial genomes (Figure 3). The results show that the order and strand orientation of the PCGs are highly conserved among genomes, with overall synteny similarity frequently reaching ≥0.95 and no large-scale rearrangements that would alter the ordering of coding regions. The differences observed are largely confined to noncoding hypervariable regions or a few tRNA clusters, manifested as changes in intergenic spacer length and minor shifts at annotation boundaries. Together with the previously noted AT bias and RSCU skew, these findings indicate that the architectural framework of the mitochondrial genome is broadly stable, whereas evolutionary changes at the sequence and codon-usage levels are more likely to drive differentiation among populations or geographic lineages. Accordingly, structural rearrangements are unlikely to be the primary mechanism underlying mitochondrial diversification in this species, and subsequent comparative and selection analyses can focus—without confounding from structural change—on nucleotide substitutions, codon-usage bias, and the ratio of nonsynonymous to synonymous substitutions (dN/dS). A comparative analysis of the lengths of the 13 mitochondrial protein-coding genes (PCGs) across 18 mitogenomes revealed highly consistent gene sizes among samples, with only minor variation detected in a few loci (notably ATP8, ND3, and ND6). ND5 was the longest gene, whereas ATP8 was the shortest, and the lengths of the COX subunits and CYTB were particularly invariant. These patterns indicate that the coding regions are under strong functional constraint and are overall conserved; the small length differences observed are likely attributable to short indels within AT-rich segments and to differences in start/stop codon usage or annotation boundaries. Consistent with our observations of AT bias and codon-usage skew, these results suggest that mitochondrial divergence in this taxon is driven more by subtle sequence and codon-level evolution than by changes in genome architecture (Figure 4).

The lrRNA and srRNA of A. cerana-Southern China were 1324 bp and 773 bp in length, respectively. lrRNA and srRNA were typically separated by trnV (Table 2). The A + T and G + C content of rRNA genes were 82.36% and 17.64%, and the AT skew and GC skew were −0.01 and −0.32, respectively, suggesting an apparent bias toward the use of T and C. Based on the length comparison (Figure 5), the srRNA is predominantly 772–787 bp and the lrRNA 1299–1337 bp, with minimal variation across the 18 mitogenomes, further supporting strong structural conservation of the rRNA region; the two rRNAs are typically separated by trnV.

In the mitochondrial genome of A. cerana (Southern China population), twenty-two tRNA genes were identified. Fourteen are encoded on the J-strand (major strand) and eight on the N-strand (minor strand). Except for trnS1, which lacks the dihydrouridine (DHU) arm, all tRNAs exhibit canonical cloverleaf secondary structures (Figures 6,7). The tRNA complement spans 1,468 bp in total, accounting for 9.10% of the mitogenome, with individual gene lengths ranging from 62 bp (trnS1) to 79 bp (trnP) (Table 3). Across genomes, tRNA lengths are tightly clustered—mostly 65–72 bp—with only minor 1-3 bp differences at a few loci, most notably in trn-Pro and trn-Thr; this pattern is consistent with small indels in AT-rich segments, subtle annotation-boundary shifts, and the missing DHU arm in trnS1 (Figure 8). These regions are strongly AT-rich (A + T = 87.35%, G + C = 12.65%; A + T is 6.91-fold higher than G + C). The AT-skew and GC-skew are 0.03 and −0.13, respectively, indicating a slight bias toward A over T and a pronounced bias toward C over G.

The AT-rich region of the A. cerana-Southern China mitochondrial genome extended over 818 bp and was located between srRNA and trnS1. The AT-rich region contained the highest A + T content (97.43%) in the entire mitogenome. Both the AT skew and GC skew for the AT-rich region were negative, indicating that T and C are more abundant than A and G, respectively (Table 3).

Phylogenetic analysis was conducted using 13 PCGs with 17 closely related taxa (Figure 9). The results of the evolutionary tree analysis revealed that A. cerana was divided into four branches, one from Taiwan and one from Borneo, with high node support. A. cerana from Korea, Japan, and Russia formed a single colony, while those from mainland China formed an additional colony. A. cerana-Southern China and other A. cerana in mainland China represent sister groups. The genetic distance of A. cerana-Southern China was closer to that of A. cerana-Central China and A. cerana-Aba, and the differentiation time for A. cerana-South Yunnan was earlier.

The circular Circos visualization shows mitochondrial homologous fragments mapping, via multiple links, to numerous loci on chr1–chr16, forming a multi-chromosomal, dispersed distribution rather than a continuous large syntenic block. The hits are predominantly short tracts (thin links with no extensive overlap), with relatively denser connections on chr1, chr6, and chr15–16. No long, one-to-one syntenic structure of mitochondrial coding regions was detected in the nuclear genome. This distribution pattern is consistent with nuclear mitochondrial DNA segments (NUMTs) (Figure 10).