High-quality Genome Assemblies of Seven Korean Honeybee Breeding Lines

1. Sample collection and preservation

Seven honeybee strains were selected from the Korean National Seed Bank, comprising five A. mellifera strains and two A. cerana strains. The A. mellifera strains included three of Italian origin (Am-A, Am-C, Am-F), one Carpatica/Carniolan-type strain (Am-D), and one Carniolan strain (Am-V). The A. cerana strains (Ac-R, Ac-X) were both of Korean origin. For each strain, more than 50 individual worker bees were collected to ensure adequate representation of colony genetic diversity. Samples were collected from two primary locations in Jeollabuk-do Province, South Korea: the Wido Isolated Breeding Station in Buan-gun and the National Institute of Agricultural Sciences in Jeonju. The Wido station is situated on an island, providing geographic isolation that minimizes the risk of genetic contamination from external colonies. All samples used for genomic analyses were nurse bees aged approximately 3-7 days post-eclosion, collected directly from the brood nest area of each colony. Immediately upon collection, specimens designated for DNA extraction were individually flash-frozen using portable liquid nitrogen dry shippers and transported to the laboratory, where they were stored at -80℃ until DNA extraction.

2. High Molecular Weight(HMW) DNA extraction

HMW DNA was extracted from whole body of individual nurse bees using a modified CTAB (cetyltrimethylammonium bromide) protocol optimized for insect tissues. Briefly, whole body were homogenized in liquid nitrogen, incubated in CTAB extraction buffer (2% CTAB, 1.4 M NaCl, 100 mM Tris-HCl pH 8.0, 20 mM EDTA, and 1% beta-mercaptoethanol) at 37℃ for 1 hour, followed by chloroform:isoamyl alcohol (24 : 1) extraction and ethanol precipitation. DNA pellets were washed with 70% ethanol, air-dried, and resuspended in TE buffer (10 mM Tris-HCl pH 8.0, 1 mM EDTA). Extract quality was assessed by 1% agarose gel electrophoresis and quantified with the Quantus Fluorometer (Promega, USA).

3. PacBio HiFi sequencing

The long read libraries were prepared using the PacBio SMRTbell Prep Kit 3.0 and SMARTbell barcoded adapter plate 3.0 according to the manufacturer’s protocol (Pacific Biosciences). Briefly, 5 μg of HMW DNA per sample was sheared to a target size of 15-20 kb using a g-TUBE (Covaris), followed by DNA damage repair, end repair, and ligation of SMRTbell adapters. Size selection was performed using the PippinHT system (Sage Science) to enrich for fragments >10 kb. Library quality and concentration were assessed using the Femto Pulse systems (Agilent) and Quantus Fluorometer. To maximize sequencing efficiency and cost-effectiveness, libraries from 7-8 individual honeybee strains were pooled in equimolar ratios and loaded onto a single SMRT Cell. Longread sequencing was performed on the PacBio Revio system using Sequencing Kit 3.0 and the 30-hour movie acquisition protocol. This pooling strategy was designed to achieve approximately 20× genome coverage per strain, assuming a genome size of approximately 230 Mbp for A. mellifera and 220 Mbp for A. cerana. Circular Consensus Sequencing (CCS) was performed using SMRT Link v13.3 software with default parameters to generate HiFi reads with a minimum predicted accuracy of Q20 (99% accuracy).

4. Oxford Nanopore Pore-C sequencing

Pore-C sequencing based on the Oxford Nanopore sequencing platform was employed to acquire chromosomal contact information for insect tissues. For each strain, three whole nurse bees were fixed with 1% formaldehyde for 10 minutes at room temperature to crosslink chromatin, followed by quenching with 125 mM glycine. Fixed cells were treated with permeabilisation solution (10 mM Tris-HCl, 10 mM NaCl, 0.2% IGEPAL CA-630). In situ digestion was performed using a cocktail of restriction enzymes (NlaIII) to fragment chromatin while preserving crosslinked contacts. Digested DNA ends were ligated using T4 DNA ligase (New England Biolabs). Following ligation, crosslinks were reversed by overnight incubation at 65℃ with proteinase K, and DNA was purified using phenol-chloroform extraction and ethanol precipitation. Purified Pore-C DNA was used to prepare Oxford Nanopore sequencing libraries using the Ligation Sequencing Kit (SQK-LSK114) according to the manufacturer’s protocol. Libraries were loaded onto Prometh ION flow cells and sequenced on a PromethION24 instrument. Sequencing was performed for 72 hours per flow cell to achieve approximately 20× genome coverage per strain. Base calling was performed using Dorado v7.9.8 in Super-accurate mode.

5. Genome assembly pipeline

Genome assembly was performed using a multi-stage pipeline integrating PacBio HiFi reads for contig assembly and ONT Pore-C data for chromosome-level scaffolding. First, HiFi reads were assembled into contigs using NextDenovo v2.5.2 with default parameters (Hu et al., 2024). NextDenovo employs a correct-then-assemble strategy for long-read sequencing data to generate highly contiguous genome assemblies. Contig-level assemblies were polished using the original HiFi reads with Next Polish v1.4.1 to correct any residual errors (Hu et al., 2020). Polished contigs were then scaffolded to chromosome level using Pore-C data. Multi-way chromatin contacts were identified using pore-c-snakemake pipeline, which filters spurious contacts, normalizes contact frequencies, and generates contact matrices. Chromosome-scale scaffolding was performed using 3D-DNA v190716, which utilizes chromatin contact frequency information to detect mis-joins and to order and orient contigs into chromosome-scale scaffolds (Dudchenko et al., 2017). Manual curation of scaffolds was performed using Juicebox Assembly Tools v2.17 to correct mis-joins and resolve ambiguities based on visual inspection of Pore-C contact maps (Dudchenko et al., 2018). The curated assembly was finalized using the 3D-DNA postreview pipeline.

6. Quality assessment

Assembly quality was assessed using multiple complementary metrics. Contiguity was evaluated using standard statistics including total assembly size, number of scaffolds, N50 (the length of the shortest scaffold in the set of scaffolds that together represent at least 50% of the assembly), N90, maximum scaffold length, and average scaffold length. These metrics were calculated using QUAST v5.2.0 (Mikheenko et al., 2018). Completeness was assessed using BUSCO (Benchmarking Universal Single-Copy Orthologs) v5.8.3 with the apoidea_odb12 database, which contains 6,507 conserved single-copy orthologs expected to be present in Hymenoptera genomes. BUSCO classifies genes as complete (single-copy or duplicated), fragmented, or missing, providing a quantitative measure of assembly completeness and potential contamination (indicated by elevated duplication rates) (Seppey et al., 2019). Pore-C contact maps were generated using 3D-DNA v190716 and visualized using Juice box v2.17. Contact maps provide a visual representation of chromatin interaction frequencies across the genome, with strong diagonal blocks indicating contiguous chromosome-scale scaffolds and off-diagonal signals revealing structural features such as centromeres, telomeres, and large-scale rearrangements.

7. Gene ontology analysis

Functional annotation of the predicted gene models for all seven chromosome-level assemblies was performed using the OmicsBox platform (v3.2.9, Biobam, Spain), integrating sequence homology searches and domain signature identifications (Götz et al., 2008). First, proteincoding gene models were predicted using AUGUSTUS (v3.5.0) implemented in the OmicsBox under ab initio prediction mode, utilizing prediction models pre-trained specifically on Apis-specific gene structures. While direct transcriptomic evidence was not integrated during the prediction step, these ab initio models were subsequently validated and filtered through extensive sequence homology searches and domain consensus analyses. The resulting coding sequences (CDSs) were subsequently queried against the NCBI non-redundant (nr) protein database using BLASTx with an e-value significance threshold of 1×10^-5 and a maximum of 20 database hits retained per sequence (Altschul et al., 1997; Sayers et al., 2025). Concurrently, functional domains, family signatures, and active sites were identified in silico by running InterProScan (v5.72-103.0) against multiple signature databases, including NCBIfam, SFLD, PANTHER, HAMAP, ProSite Profiles, CDD, PRINTS, Pfam, PIRSF, SUPERFAMILY, Gene3D, AntiFam, FunFam, and PIRSR (Jones et al., 2014). Gene Ontology (GO) terms representing three primary domains-Biological Process (BP), Molecular Function (MF), and Cellular Component (CC)-were assigned by merging the BLASTp mapping results and InterProScan domain annotations within the OmicsBox suite (Consortium, 2004). Enzyme Commission (EC) numbers were also mapped to sequences containing valid metabolic enzyme signatures. To perform high-resolution functional profiling, direct GO term counts (capturing levels 3 to 7+) were extracted from the annotation outputs for each strain (e.g., go_count_bp.txt, go_count_mf.txt, and go_count_cc.txt) to examine granular functional differences. For comparative genomic analyses, GO level 2 distribution counts were normalized as percentages of the total GO-annotated genes within each strain to account for variations in gene-space prediction sizes. To determine if the overall functional landscape varied significantly among the seven honeybee breeding lines, a Chi-square (χ²) test of independence was applied to the GO term contingency table. Subsequently, a two-tailed Welch’s t-test (independent t-test with unequal variances) was performed to identify specific GO categories exhibiting statistically significant differences (p<0.05) between A. cerana (n=2) and A. mellifera (n=5). Finally, a Z-score outlier analysis was conducted across the seven strains to identify unique, lineage-specific functional expansions (Z>2.0). All statistical calculations, Welch’s t-tests, and high-resolution multi-panel bar charts and heatmaps were executed in Python (v3.10) using scipy.stats, matplotlib, and seaborn libraries.

8. Phylogenetic analysis

To confirm the genetic identity and evolutionary relationships of the seven sequenced strains, we performed phylogenetic analysis based on complete mitochondrial genomes. The consensus contigs were assembled by NextDenovo (v2.5.2) from the selected HiFi reads corresponding to the bee reference mitochondrial genome, followed by concatenation using in-house perl script. The draft mitocondrial genome was completed to the circular form and polished using HiFi long reads by NextPolish (v1.4.1). The assembled mitochondrial genome was set to start from the same location with NC_014295.1 (A. cerana) or NC_051932.1 (A. mellifera). Mitochondrial genome sequences from the seven Korean strains were aligned with publicly available mitogenomes representing major Apis species and A. mellifera subspecies, including A. m. ligustica, A. m. carnica, A. m. mellifera, A. m. anatoliaca, and A. cerana japonica. To provide an evolutionary context, we included several Apis species as reference taxa (A. florea, A. andreniformis, A. dorsata, A. laboriosa). Multiple sequence alignment was performed using MAFFT v7.526 with the FFT-NS-2 algorithm. Phylogenetic trees were inferred using maximum likelihood (ML) in IQ-TREE v2.4.0 with 1,000 ultrafast bootstrap replicates to assess node support (Nguyen et al., 2015). The best-fit substitution model was selected based on the Bayesian Information Criterion (BIC). Trees were rooted using Bombus ignitus (GenBank accession: DQ870926) as the outgroup and visualized using FigTree v1.4.5 and MEGA11 (Katoh and Standley, 2013).

1. Geographic isolation and genetic background of conserved Apis breeding lines

Seven honeybee strains were successfully collected from two primary locations in Jeollabuk-do Province, South Korea: the Wido Isolated Breeding Station in Buan-gun and the National Institute of Agricultural Sciences in Jeonju. The five A. mellifera strains represented three major European subspecies: Italian bees (A. m. ligustica: Am-A, Am-C, Am-F), Carpatica/Carniolan-type bees (Am-D), and Carniolan bees (A. m. carnica: Am-V). The two A. cerana strains (Ac-R, Ac-X) were native Korean honeybees. All strains are registered in the Korean National Seed Bank and maintained under controlled breeding conditions to preserve genetic purity and prevent introgression from external populations. More than 50 individuals per strain were collected to ensure adequate sampling of within-colony genetic diversity.

2. High-depth long-read sequencing metrics and comparative performance of de novo genomic assemblers

PacBio HiFi sequencing on the Revio platform generated high-quality long-read data for all seven strains (Supplementary Table 1). The total number of HiFi reads per strain ranged from 777,008 (Am-C) to 1,116,667 (Am-V), corresponding to total yields of 10.35-12.45 Gbp. Read length N50 values ranged from 10,247 bp (Ac-R) to 17,038 bp (Am-C), with a median of approximately 13-14 kb for most strains. Median read quality scores were exceptionally high, ranging from Q36 to Q39, with 100% of reads exceeding Q20 (99% accuracy), 69-80% exceeding Q30 (99.9% accuracy), and 28-46% exceeding Q40 (99.99% accuracy). These quality metrics confirm the suitability of the HiFi data for high-fidelity genome assembly. The sequencing depth for each strain was estimated at approximately 45-55× coverage based on total yield and expected genome sizes (approximately 230 Mbp for A. mellifera and 220 Mbp for A. cerana). This exceeds the target coverage of 20× and provides substantial redundancy for accurate contig assembly and variant calling.

To evaluate assembly quality across the seven bee strains, we performed de novo assembly using two different long-read assemblers: NextDenovo and Hifiasm (Table 1). The resulting contig-level statistics showed distinct differences in total assembly size, contig numbers, and N50 lengths between the two methods. For the A. mellifera strains, NextDenovo generated assemblies with total sizes ranging from 226.96 Mb (Am-F) to 233.15 Mb (Am-D), comprised of 65 to 91 contigs. The contig N50 values for these assemblies ranged from 9.03 Mb (Am-D) to 11.63 Mb (Am-A). In comparison, Hifiasm produced larger assembly sizes for the same strains, ranging from 264.31 Mb (Am-F) to 278.97 Mb (Am-A), with higher contig counts ranging from 287 to 486, although it yielded contig N50 values between 12.31 Mb (Am-C) and 14.23 Mb (Am-A). In the A. cerana strains, the variance between the two assemblers was more pronounced. NextDenovo assembled the Ac-R and Ac-X genomes into 130 contigs (Size: 216.05 Mb, N50: 6.00 Mb) and 62 contigs (Size: 218.37 Mb, N50: 10.63 Mb), respectively. Conversely, Hifiasm assemblies for Ac-R and Ac-X resulted in substantially larger total sizes of 513.74 Mb and 508.36 Mb, with contig numbers increasing to 4,207 and 1,875, and contig N50 values of 1.21 Mb and 1.78 Mb, respectively. Based on the total assembly sizes and contig continuity metrics, the NextDenovo assemblies were selected as the representative genomes for all seven strains and utilized for subsequent downstream analyses.

Table 1.

Comparative summary of genome assembly metrics between NextDenovo and Hifiasm pipelines

3. Pore-C facilitated chromosomal scaffolding and species-specific scaffold length polymorphisms

The final chromosome-level assemblies captured a total length ranging from 207.68 Mb (Ac-R) to 223.22 Mb (Am-C) (Table 3). For all seven strains, the genomic sequences were successfully anchored into 16 major pseudomolecules, which precisely corresponds to the haploid chromosome number (n=16) of the genus Apis (Fig. 1, Supplementary Table 2). The genome-wide Pore-C interaction heatmaps demonstrated highly organized diagonal contact intensities with 16 distinct interaction blocks, validating the structural accuracy and continuity of the chromosomal scaffolding. The length distribution of the 16 chromosome-level scaffolds exhibited high consistency within each species, while displaying distinct speciesspecific profiles between A. mellifera and A. cerana. The five A. mellifera strains (Am-A, Am-C, Am-D, Am-F, and Am-V) showed nearly identical chromosomal length profiles across all lineages. In these strains, Scaffold_1 represented the largest chromosome, with lengths consistently clustered between 27.75 Mb and 27.79 Mb, whereas the smallest chromosome, Scaffold_16, ranged from 7.58 Mb to 7.63 Mb. The intermediate pseudo-molecules maintained highly uniform sizes among the five lineages, indicating strong structural conservation within the evaluated A. mellifera genomes. In contrast, the two A. cerana strains (Ac-R and Ac-X) presented intra-species variation in specific scaffold sizes compared to A. mellifera. In the Ac-R genome, Scaffold_1 measured 21.05 Mb, which was notably shorter than the 26.97 Mb observed in Ac-X. Conversely, Scaffold_2 and Scaffold_3 in Ac-R were 20.59 Mb and 18.01 Mb, respectively, both of which were larger than the corresponding scaffolds in Ac-X, which measured 17.17 Mb and 15.84 Mb. For the remaining pseudo-molecules from Scaffold_4 through Scaffold_16, both strains maintained comparable length distributions, with Scaffold_16 terminating at approximately 7.05 Mb for Ac-R and 7.26 Mb for Ac-X. These chromosomal statistics and the clear separation of the 16 contact matrices confirm that the assemblies achieved high-fidelity, chromosome-scale resolution suitable for subsequent downstream investigations.

Fig. 1.

Genome-wide Pore-C interaction heatmaps of seven reconstructed honeybee (Apis) strains. The individual contact matrices illustrate the chromatin interaction densities across the chromosome-level genome assemblies for five A. mellifera lineages-Am-A, Am-C, Am-D, Am-F, and Am-V-and two A. cerana lineages-Ac-R and Ac-X. The color intensity or pixel density reflects the frequency of localized physical genomic interactions, with the brightest diagonal signals representing contiguous genomic blocks. For all seven strains, the contact matrices display 16 sharply resolved, discrete interaction blocks along the diagonal axes. This configuration aligns precisely with the haploid chromosome number (n=16) characteristic of the genus Apis, validating the structural completeness, high continuity, and robust anchoring of the de novo scaffolded pseudo-molecules without significant inter-chromosomal misassembles.

4. Assessment of biological completeness and gene-space quality via ortholog verification

To evaluate the biological completeness and gene-space quality of the final scaffolded assemblies, a BUSCO analysis was performed using a core set of 6,507 orthologs (Table 2). The assessment demonstrated exceptionally high and uniform completeness across all seven evaluated bee genomes, with no significant structural or genecontent variation observed between the species or strains. For the A. mellifera strains, the total percentage of complete BUSCOs was remarkably consistent, ranging from 98.1% in Am-A, Am-C, Am-D, and Am-F to 98.2% in Am-V. Within these complete orthologs, single-copy genes accounted for 98.0% to 98.1% of the total dataset, while duplicated BUSCOs remained extremely low at a mere 0.1% (4 to 6 genes) for all five strains. The fraction of fragmented and missing BUSCOs was also tightly constrained, spanning only 0.5% (30 to 31 genes) and 1.4% to 1.5% (89 to 95 genes), respectively. A similarly high level of genomic completeness was observed in the A. cerana strains, with both Ac-R and Ac-X achieving a total complete BUSCO rate of 98.1%. Single-copy BUSCOs reached 98.0% for Ac-R and 97.9% for Ac-X, accompanied by a low duplication rate of 0.1% (7 and 8 genes, respectively). The fragmented genes constituted 0.5% to 0.6%, and the missing orthologs accounted for a minor fraction of 1.4% in both genomes.

Table 2.

BUSCO assessment results of genome assemblies for Apis mellifera and Apis cerana strains using the apoidea_odb12 lineage dataset

5. Structural gene prediction fidelity and functional Gene Ontology (GO) profiling

Following the chromosome-level scaffolding, structural gene prediction and functional annotation were performed to evaluate the coding potential of the seven assembled genomes (Table 3). The total number of predicted protein-coding genes ranged from 6,860 (Ac-X) to 6,915 (Ac-R) in A. cerana, and from 7,261 (Am-F) to 8,751 (Am-D) in A. mellifera. To evaluate the accuracy of the de novo gene prediction models, the predicted gene sets were cross-referenced against multiple public functional databases, exhibiting exceptionally high homology rates. Across all seven strains, approximately 97% to 98% of the total predicted genes were successfully validated through alignment against the NCBI non-redundant (nr) database. Specifically, the nr mapping rates spanned from 97.20% (Ac-X) to 97.24% (Ac-R) for A. cerana, and from 96.94% (Am-A) to 97.58% (Am-V) for A. mellifera. Furthermore, protein domain analysis using InterProScan (IPS) successfully identified known functional domains in 65.71% to 66.22% of the A. cerana genes and 66.56% to 67.64% of the A. mellifera genes. The integrated functional annotation pipeline culminated in successful Gene Ontology (GO) term assignment. Through the merging of IPS and blast-derived annotations, approximately 90% of the entire predicted gene complements were successfully assigned definitive GO terms. The merged GO annotation rates were highly uniform across all lineages, ranging from 89.75% to 90.45% in A. cerana and from 89.40% to 90.13% in A. mellifera. This high proportion of functionally characterized genes, combined with the comprehensive database alignment rates, demonstrates that the de novo gene prediction pipeline achieved substantial structural integrity without major sequence fragmentation or gene loss, confirming the high quality of the final genome annotations.

Table 3.

Summary of genome functional annotation metrics for Apis mellifera and Apis cerana strains

To further characterize the functional distribution of the annotated genes, the assigned GO terms were classified into Level 2 categories under three main ontologies: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) (Fig. 2). The overall distribution of GO terms exhibited a highly synchronized and consistent pattern across all seven evaluated strains. Within the Biological Process category, cellular process emerged as the most predominant term, encompassing approximately 76% of the GO-annotated genes in all lineages. This was followed by biological regulation and regulation of biological process, which consistently accounted for roughly 27-29% and 26-28% of the functional repertoire, respectively. Other well-represented BP subcategories included localization (~18%) and response to stimulus (~15%), demonstrating a highly uniform allocation of biological roles across both A. mellifera and A. cerana genomes. In the Molecular Function domain, binding and catalytic activity constituted the most dominant functional assignments. The binding category represented the largest fraction, with a stable proportion of approximately 83% across all seven strains. Catalytic activity was the second most abundant function, ranging between 39% and 47% of the annotated coding space. The remaining MF subcategories, such as transporter activity, transcription regulator activity, and ATP-dependent activity, presented identical hierarchical rankings and minor percentage variations among the lineages. For the Cellular Component ontology, the annotated genes were heavily concentrated in two primary structural categories: cellular anatomical structure, which represented over 80% of the annotated pool, and protein-containing complex, which comprised approximately 21-23% of the genes across all strains. The negligible divergence in GO term densities and distribution patterns across the entire dataset underlines that despite the genomic variance and lineage-specific divergence observed in structural statistics, the fundamental core proteome and functional machinery remain strictly conserved between A. mellifera and A. cerana.

Fig. 2.

Gene Ontology (GO) Level 2 functional classification of annotated protein-coding genes across seven Apis strains. The bar charts illustrate the distribution and percentage of functional categories annotated within the genomes of five A. mellifera strains (Am-A, Am-C, Am-D, Am-F, and Am-V) and two A. cerana strains (Ac-R and Ac-X). The classification is divided into three primary ontologies: (A) Biological Process, (B) Molecular Function, and (C) Cellular Component. The horizontal axis indicates the percentage of total GO-annotated genes within each individual lineage. Color-coded bars represent specific honeybee lineages, as defined in the species legend box. Comparative distribution displays a highly conserved and uniform functional landscape across all evaluated genomes, highlighting the stable retention of core physiological and metabolic machinery within the genus Apis.

6. Comparative functional enrichment and strain-specific toxicological and adaptive outliers

To capture functional variations within the conserved genomic architecture, a comparative enrichment analysis of GO Level 2 terms was conducted using Fisher’s exact test, followed by Benjamini-Hochberg False Discovery Rate (FDR) correction (adjusted p<0.05), based on the pooled counts of annotated genes across all genomes (Fig. 3A, Supplementary Table 3). A. cerana genomes exhibited significant enrichment in pathways associated with cellular anatomical structure (+2.88%, adjusted p<0.001) and protein-containing complexes (+2.02%, adjusted p<0.001). This structural profile was further associated with enrichment in biological regulation (+1.57%, adjusted p=0.005), developmental processes (+0.84%, adjusted p=0.007), and regulation of biological processes (+1.45%, adjusted p=0.009). Conversely, A. mellifera genomes were significantly enriched in categories driving metabolic efficiency and transport. The most pronounced divergence was observed in catalytic activity, which showed a significant difference of 5.42% relative to A. cerana (adjusted p<0.001). Additionally, transporter activity was significantly enriched in A. mellifera (2.10%, adjusted p<0.001), along with ATP-dependent activity (1.18%, adjusted p<0.001) and localization (1.10%, adjusted p=0.024).

Fig. 3.

Comparative enrichment and strain-specific outlier analyses of GO Level 2 terms among Apis lineages. (A) Species-specific functional enrichment under Fisher’s exact test (FDR-adjusted p<0.05). The horizontal axis displays the difference in the relative percentage of annotated genes between A. cerana and A. mellifera (MeanAc - MeanAm). Red bars extending to the right represent functional categories significantly enriched in A. cerana, while blue bars extending to the left indicate categories significantly overrepresented in A. mellifera. Statistical significance levels (p-values) and mean difference values are explicitly labeled on individual bars. (B) Outlier functional adaptations of individual strains based on Z-score profiling across the seven genomes. The horizontal axis represents the standardized deviation (Z-score) of GO term proportions for each strain against the group mean. Vertical red dashed lines denote the statistical outlier thresholds (|Z|=2.0). Color-coded bars mirror the specific honeybee lineages indicated in the color key. Prominent lineage-specific signatures are highlighted, including the acute elevation of toxin activity in the Am-V strain (Z=2.45), revealing highly localized and specialized evolutionary selection pressures acting upon individual genomic landscapes.

Beyond the species-level evolutionary divergence, a standardized Z-score analysis (|Z|>2.0) was deployed to screen for strain-specific functional outliers across the seven evaluated lineages (Fig. 3B). This filtration identified localized genomic configurations unique to individual bee strains that deviated significantly from the crosslineage baseline. A distinct functional outlier signature was captured in the A. mellifera Carniolan-derived strain, Am-V, which displayed a pronounced statistical overrepresentation specifically within the toxin activity category (Z=2.45). Furthermore, the Z-score profiling revealed positive deviations shared between the two Carniolan lineages, Am-D and Am-V, in specific somatic and metabolic categories. In the Z-score distribution of these two strains, this profile was characterized by the statistically significant overrepresentation of crucial metabolic and xenobiotic detoxification gene families, including cytochrome P450 monooxygenases. The quantitative convergence of these amplified detoxification pathways and metabolic homeostasis categories was uniquely localized to the Am-D and Am-V datasets compared to the other three western honeybee lineages.

7. Mitochondrial genome characterization and phylogenetic placement

To determine the maternal lineages and evolutionary origins of the seven honeybee strains, we successfully reconstructed their complete, gapless circular mitochondrial genomes utilizing the high-fidelity long reads from the PacBio Revio platform. The assembled mitogenomes exhibited highly conserved structural architectures and species-specific length distributions (Supplementary Table 4). The two A. cerana strains (Ac-R and Ac-X) possessed concise genomes ranging from 16,224 bp to 16,247 bp. In comparison, the five A. mellifera strains (Am-A, Am-C, Am-D, Am-F, and Am-V) displayed slightly larger but exceptionally uniform sizes, spanning from 16,515 bp to 16,576 bp. All seven mitogenomes exhibited an extreme AT-bias, with the total AT content ranging between 84.29% and 85.09%.

To evaluate the phylogenetic positions of these lineages, a comprehensive 65-taxon Maximum Likelihood (ML) phylogenetic tree was constructed based on the complete mitogenome sequences (Fig. 4, Supplementary Table 5). The phylogenetic topology revealed a clear genetic dichotomy that separated the eastern and western honeybee lineages into distinct, well-supported geographic clades. The two Korean A. cerana genebank strains (Ac-R and Ac-X) clustered together with absolute bootstrap support within the mainland East Asian continental clade (y=1.0 to 9.0 region), positioning adjacent to reference accessions from China and India. Notably, this mainland continental clade was phylogenetically distinct from the Japanese A. cerana lineage (AP017985), as well as from other wild genetic resources native to Southeast and Southwest Asia. Conversely, the five A. mellifera strains maintained by the national genebank (Am-A, Am-C, Am-D, Am-F, and Am-V) were unequivocally assigned to the European C-lineage (y=23.0 to 38.0 region), which encompasses major Mediterranean and Eastern European subspecies such as A. m. carnica and A. m. ligustica. Although all five strains anchored firmly within the C-lineage, they exhibited subtle branch-length variations and distinct sub-clade bifurcations, revealing a meaningful degree of intra-lineage genetic diversity. Particularly, the Am-C strain clustered closely with Eastern Mediterranean gene pools, including the A. m. anatoliaca subspecies, in addition to the typical Western European references.

Fig. 4.

Reconstructed global phylogenomic tree of the genus Apis and its subspecies based on complete mitogenomes. The Maximum Likelihood (ML) phylogenetic topology was inferred using IQ-TREE 2 based on the alignment of 65 complete mitochondrial genomes via MAFFT. Bombus ignitus (GenBank accession: DQ870926) was utilized as the evolutionary outgroup to root the tree. The horizontal scale bar at the bottom indicates the evolutionary distance, measured in substitutions per site. The red shaded box highlights the A. cerana continental clade (y=1.0 to 9.0 region), positioning the two Korean seed bank strains (Ac-R and Ac-X, marked with asterisks) strictly within the mainland East Asian lineage alongside Chinese and Indian accessions. The blue shaded box outlines the A. mellifera evolutionary lineage (y=23.0 to 38.0 region), mapping the five Korean seed bank strains (Am-A, Am-C, Am-D, Am-F, and Am-V, marked with asterisks) exclusively into the European C-lineage group.

To overcome the maternal-inheritance limitation of mitochondrial phylogeny and confirm nuclear breeding-line integrity, we also performed a whole-genome nuclear distance comparison using Mash (k=21, s=10,000) (Supplementary Fig. 1). The nuclear distance matrix revealed a clear genetic divergence between A. cerana and A. mellifera strains (average distance of 0.058, equivalent to 94.2% Average Nucleotide Identity [ANI]). Within species, the nuclear pairwise distances were tightly constrained (~0.005 distance, 99.5% ANI), with the two Korean A. cerana seed bank strains (Ac-R and Ac-X) clustering closely together. This dual phylogenomic validation confirms the genetic purity, breed integrity, and maternal-paternal lineage stability of all seven national seed bank repositories.

1. Chromosome-scale length variations and implications for population divergence

A compelling structural finding in our genomic survey was the intra-species scaffold length variance observed between the two A. cerana lineages, Ac-R and Ac-X. Specifically, scaffold_1 in Ac-R was found to be approximately 5.92 Mb shorter than its homologue in Ac-X, whereas scaffold_2 and scaffold_3 in Ac-R exhibited reciprocal size increases of 3.42 Mb and 2.16 Mb, respectively. It is critical to underscore that while formal structural variant (SV) discovery pipelines-such as long-read split-mapping or optical mapping-were not deployed in the current study, these macro-level length disparities captured via Pore-C scaffolding provide indirect yet tangible genomic signatures of large-scale chromosomal rearrangements. In hymenopteran genomics, such pronounced chromosomal length discrepancies among closely related allopatric or sympatric populations typically point to massive structural modifications, most notably chromosomal inversions, translocations, or the localized lineage-specific accumulation of transposable elements (TEs) (Kapusta et al., 2017; Mérot et al., 2020). Large chromosomal inversions, for instance, frequently act as evolutionary “supergenes” by physically suppressing meiotic recombination between divergent lineages (Wellenreuther and Bernatchez, 2018). This recombination barrier locks adaptive allele complexes together, allowing populations to maintain distinct ecological phenotypes-such as specialized overwintering strategies, foraging behaviors, or predatory defense mechanisms-despite potential gene flow (Wallberg et al., 2017).

Given that Ac-R and Ac-X represent distinct indigenous conservation lines within the Korean peninsula, the identified polymorphic regions on scaffolds 1, 2, and 3 stand out as primary candidates for maintaining unique agronomic traits. To transition from these macro-structural observations to definitive mechanistic insights, a comprehensive, high-resolution SV profiling initiative utilizing long-read graph-pangenome mapping remains an essential next step (Koren et al., 2018; Garg et al., 2021). Characterizing the exact breakpoints and sequence content of these polymorphic blocks will be paramount to unlocking the genetic basis of fine-scale adaptive divergence in Korean A. cerana, directly empowering marker-assisted precision breeding and regional ecotype preservation (Meuwissen et al., 2001; Brascamp et al., 2018).

2. Species- and strain-specific functional divergence and phenotypic correlates

The comparative statistical framework identified targeted functional shifts that suggest potential genomic underpinnings for long-observed phenotypic, behavioral, and agronomic disparities between eastern and western honeybees, although direct links between broad GO Level 2 terms and complex phenotypes require cautious interpretation. The enrichment of biological regulation and regulation of biological process in A. cerana may reflect its complex regulatory adaptations. Historically, A. cerana has demonstrated remarkable ecological resilience against devastating endemic parasites and pathogens, such as Varroa mites and the Sacbrood virus, largely driven by highly specialized hygienic behaviors including intensive grooming and targeted removal of infected brood (Boecking and Genersch, 2008; Oxley et al., 2010). While broad GO categories such as response to stimulus or immune system process do not show significant species-wide expansion in overall gene counts under strict FDR correction, the specific gene families nested within these categories-including chemosensory receptors and specific immune signaling cascades detailed in subsequent sections-provide the precise molecular repertoire supporting these rapid, systemic responses. Furthermore, the robust representation of cellular and regulatory pathways provides a plausible genetic background that may facilitate the acute behavioral vigilance of A. cerana, potentially modulating rapid colony-level alarm signaling and collective defensive tactics-such as heat-balling behaviors-when encountering aggressive native predators like Vespa mandarinia (Tan et al., 2012).

Conversely, the genomic architecture of A. mellifera appears evolutionarily tailored toward maximized metabolic output and resource acquisition. The pronounced enrichment in catalytic activity (adjusted p<0.001, -5.42%) and transporter activity (adjusted p<0.001, -2.10%) provides a compelling genetic framework that may underlie the high foraging efficiency and commercial productivity that define the western honeybee (Wallberg et al., 2014; Uzunov et al., 2017). These metabolic expansions suggest a potentially intensified enzymatic machinery dedicated to rapid carbohydrate processing, ATP-driven flight-energy generation, and systemic nutrient transport (Foret et al., 2012; Kapheim et al., 2015). This amplified metabolic throughput directly translates into the extended foraging ranges, prolonged flight capacities, and superior honey-storing behaviors observed in agricultural settings, reflecting the plausible historical selection pressures that favored high-yield traits in A. mellifera (Harpur et al., 2012).

At the intra-species level, the standardized Z-score filtering successfully exposed localized, hyper-specialized adaptations within individual commercial lines, most notably in the Carniolan-derived genotypes. The sharp statistical deviation in the toxin activity category observed exclusively in the Am-V strain (Z=2.45) marks it as an evolutionary outlier. Detailed sequence analysis of the genes driving this category in Am-V revealed a combination of host defensive peptides and gut symbiont sequences. Specifically, we identified omega-conotoxin-like protein 1 (also known as inhibitor cysteine knot peptide), a eukaryotic disulfide-rich peptide known to be involved in arthropod defense and venom. However, the outlier status was also driven by the presence of two co-assembled hypothetical proteins from the honeybee gut symbiont Snodgrassella alvi (Supplementary Table 6). This indicates that the apparent statistical enrichment of toxin activity in Am-V is a composite signature, reflecting both host defense genes and co-assembled symbiotic microbiome elements rather than a massive host-specific duplication of venom components. Furthermore, the distinct functional configurations shared by the Carniolan sister strains, Am-D and Am-V, provide a molecular link to their renowned overwintering capacity and cold hardiness (Büchler et al., 2014). Cross-strain copy number comparisons of the Cytochrome P450 (CYP) gene family (Supplementary Table 7) revealed highly conserved profiles across all seven strains, with 18 genes in A. cerana and 19 to 21 genes in A. mellifera (Am-A: 19, Am-C: 21, Am-D: 20, Am-F: 21, Am-V: 21), indicating that no gene family expansion occurred in the Carniolan strains. Thus, their cold-adaptation capabilities are likely regulated through transcriptional modulation of these conserved detoxifying pathways and metabolic homeostatic systems, rather than genomic copy number expansion. This conserved enzymatic buffer, combined with lineage-specific expression dynamics, equips these specific lineages to effectively mitigate low-temperature stress and cellular toxicity (Guarna et al., 2017), contributing to their survival rates during harsh wintering periods.

3. Phylogeographic origins and genetic identity of national repositories

The gapless reconstruction of complete circular mitogenomes and the subsequent 65-taxon phylogenomic mapping provide unequivocal maternal evolutionary contexts for the honeybee genetic resources conserved in South Korea. For the indigenous eastern honeybee strains, Ac-R and Ac-X, their robust positioning within the mainland East Asian continental clade (y=1.0 to 9.0) settles a critical question regarding their phylogeographic ancestry (Radloff et al., 2010). By exhibiting a clear genetic divergence from the Japanese A. cerana lineage (AP 017985) and wild island ecotypes of Southeast Asia, our findings demonstrate that the Korean native bee repository shares a direct evolutionary trajectory with the mainland Eurasian continental stock (Takahashi et al., 2016; Nunn et al., 2022). This biogeographic assignment is highly congruent with multi-locus nuclear screening demonstrating that South Korean A. cerana stocks form a highly integrated, uniform temperate cluster alongside Russian Far East populations, sharply distinguished from tropical Southeast Asian ecotypes (Ilyasov et al., 2022; Kaskinova et al., 2022). Critically, this regional gene pool bears the molecular signatures of severe evolutionary bottlenecks and intense purifying selection driven by historical Sacbrood virus epidemics that previously eliminated up to 95% of the domestic colonies, validating our conserved seed bank strains as indispensable, highly screened genetic reservoirs for regional ecotype preservation. This continental genomic signature implies a shared ancestral adaptation to severe temperate seasonal fluctuations, serving as a highly valuable genetic baseline for the preservation of pure domestic ecotypes and the development of robust, climate-resilient local breeding lines (Chen et al., 2016).

Parallel to this, the phylogenomic clustering of the five western honeybee strains (Am-A, Am-C, Am-D, Am-F, and Am-V) exclusively within the European C-lineage (y=23.0 to 38.0) confirms the structural retention of their Mediterranean and Eastern European maternal ancestry (De la Rúa et al., 2009; Boardman et al., 2019). More importantly, the identification of micro-clade bifurcations and subtle branch-length variations among these five lines reveals that they do not constitute a homogenous, bottlenecked population. The close phylogenetic affinity of the Am-C strain toward Eastern Mediterranean gene pools, such as A. m. anatoliaca, alongside typical A. m. carnica or A. m. ligustica frameworks, serves as a molecular footprint of historical genetic admixture. Rather than representing isolated, singular western imports, these national conservation lines function as dynamic genetic reservoirs (Harpur et al., 2012, 2015). Over decades of localized maintenance, artificial selection, and selective breeding within the Korean peninsula, these strains have integrated diverse elite European genetic components. This hybrid-driven genomic plasticity likely accelerated their regional acclimatization, resulting in a unique, domestic-adapted gene pool tailored to the specific environmental and apicultural demands of South Korea.