Introduction
The honey bee Apis mellifera L. plays a crucial role not only in natural ecosystems but also in the economy of Ukraine. Its importance is particularly evident in agriculture, where bees provide pollination of major crops and produce valuable products, including honey, bee bread, propolis, wax, royal jelly, and drone brood homogenate (Volkov et al., 2024). Ongoing global changes in the biosphere frequently result in adverse effects on the survival and stability of Apis mellifera populations (Kükrer et al., 2021; Toby et al., 2025; Barroso et al., 2025). This has led to the development of international initiatives aimed at maintaining sustainable bee populations (SDGs, Goal 2; Goal 15). Behavioral adaptation (Haider et al., 2025) and selective breeding (Uzunov et al., 2017; Brascamp et al., 2024) are among the key mechanisms mitigating anthropogenic pressures. In parallel, considerable research efforts are directed toward protecting bees and other insects from pathogens (Ferrufino et al., 2024; Zhang et al., 2025a; Darren et al., 2025). Morphometric traits of honey bees are also actively studied, as they are closely linked to environmental adaptation and functional performance (Zhang et al., 2025b). One of the most critical challenges remains the control of the mite Varroa destructor, which is responsible for large-scale colony losses across multiple continents (De la Mora et al., 2025; Hilsmann et al., 2025).
The origin, distribution, and phylogeny of Apis mellifera across Eurasia are generally well understood (Meixner et al., 2025). It is widely accepted that the evolutionary lineage C expanded into Europe from the Middle East. According to the classification of Ruttner (1988), A. m. macedonica, a representative of lineage C, is naturally distributed in North Macedonia, Bulgaria, and parts of Ukraine. However, recent studies indicate that this subspecies is being progressively displaced from its historical range due to the expansion of A. m. carnica (Tanasković et al., 2022). Notably, A. m. macedonica is scarcely addressed in some studies on the distribution of Apis mellifera (Han et al., 2012), likely due to the limited availability of genetic data. According to studies conducted by Ukrainian researchers, honey bee populations in the Carpathian region are attributed to A. m. carnica (Papp et al., 2015; 2017a; 2017b), whereas the rest of Ukraine is predominantly inhabited by Ukrainian steppe bees, presumably related to A. m. macedonica (Polishchuk, 2007; Metlytska et al., 2012). Nevertheless, the subspecies affiliation and taxonomic status of Ukrainian steppe bees remain unresolved (Engel, 1999).
Drones of Apis mellifera play a key role as carriers of the haploid genome, ensuring the transmission of genetic traits to subsequent generations. Consequently, numerous studies have focused on their behavioral (Neubauer et al., 2023; Hagan et al., 2024), physiological (Divasón et al., 2023; Yániz et al., 2024), and biological characteristics (Vasquez-Valencia et al., 2024). Genetic studies of A. m. caucasica populations in protected areas of Türkiye revealed low genetic diversity and high levels of inbreeding (Yıldız et al., 2023), which has important implications for breeding and conservation strategies. Morphometric characteristics of drones have also attracted increasing research attention (Neto et al., 2022; Zhang et al., 2024; Gilchrist et al., 2024). Particular emphasis has been placed on the analysis of drone wing morphology, as it provides highly informative traits for taxonomic and population-level studies compared to other morphological features (Krtinic et al., 2016). However, such studies are constrained by the seasonal presence of drones in colonies, typically limited to the summer period, which restricts sampling opportunities (Henriques et al., 2020; Tofilski et al., 2024). Recent investigations (Yarovets et al., 2023, 2024; Galatiuk et al., 2024a, b; Babenko et al., 2024a, b) have partially expanded the available data on drone wing morphometry in Ukraine. The present study builds upon these findings and aims to extend the application of drone wing morphometry for assessing interactions among subspecies, populations, ecotypes, and lines, as well as for estimating the degree of hybridization.
It is acknowledged among Ukrainian beekeepers that A. m. caucasica may be present in local apiaries, a hypothesis partially supported by scientific evidence (Meixner et al., 2007). This subspecies is widely distributed across Europe and Asia and has been introduced to other continents (Janashia et al., 2025), which necessitates further investigation of its phenotypic variability and interactions with other subspecies under different environmental conditions (Tavakoli et al., 2025). In this context, the present study also evaluates the extent of the influence of Caucasian bees on the morphological traits and genetic structure of honey bee populations in Ukraine.
The results of the morphometric analysis of drone wings are sufficiently robust and are summarized in the form of reference morphometric templates for local honey bee populations, indicating that the objectives of the study have been successfully achieved.
Materials and Methods
A total of 2,694 drone wings were examined from an apiary located in the settlement of Kuzemyn (on the border of the Sumy and Poltava regions), as well as 988 wings from an apiary in the settlement of Pochyno-Sofiivka (Dnipropetrovsk region).
The data processing algorithm included the following stages:
(1) digitization of wings; (2) primary classification of the entire dataset with cluster identification; (3) classification of wings from individual colonies into two clusters followed by their identification; (4) exclusion of clusters that could not be reliably identified from the dataset; (5) reclassification of the corrected dataset with the formation of morphometric templates (reference patterns) for the studied apiaries according to the methodology described in Galatiuk et al. (2024b); (6) additional reclassification of two wing clusters for each of the studied colonies (a four-cluster classifier was applied). Wing digitization involved determining the coordinates of vein intersection points according to the scheme shown in Figure 1, using the TpsDig2 software. The obtained coordinates were used to calculate morphometric traits (indices): Ci, Dbi, Pci, Disc. sh., and Ri according to Boudga et al. (2011), as well as four additional indices proposed by the authors. The last four indices were calculated as ratios of distances between landmarks on the wing: Сі.1.1=d(6;1)/d(5;4), Сі.2.1=d(10;4)/d(5;6), Сі.2.2=d(10;4)/d(8;6), Сі.3=d(9;17)/d(9;8), where d(m;n) denotes the distance between landmarks nubered m and n (Fig. 1).

Scheme of numbering and order of applying reference points on the wing.
The classification of feature arrays was performed using discriminant analysis implemented in the STATISTICA software package, with the optimal number of clusters determined according to the criterion of maximum classification correctness. As a result, a classifier was obtained and subsequently used to calculate the Degree of Similarity (DS) of cluster phenotypes to other, so-called “introgressive” subspecies and populations. The obtained DS values allow for the assessment of the Degree of Hybridization (DH) of haploid genomes of drones and queens in the component associated with wing morphology (Babenko et al., 2024a).
Subspecies and population identification of the clusters was carried out by comparison with morphometric templates obtained in previous studies (Galatiuk et al., 2023; 2024b), based on Mahalanobis distance (MD) values between the centroids of clusters and templates. The interpretation of MD values was defined as follows: 0–2 — high similarity; 2.0–2.6 — moderate similarity; 2.6–3.5 — low similarity; >3.5 — no similarity (Yarovets et al., 2024). A similar approach was applied for the identification of two wing clusters within each individual colony.
The identification of phenotypes for 22 wing clusters across 11 colonies from an apiary in the Dnipropetrovsk region was performed using the conventional approach by comparison with available morphometric templates. The assessment of the degree of hybridization of the main population with “introgressive” components was carried out using a four-cluster classifier obtained from variance analysis of data for 2,420 wings.
Results
Using discriminant analysis, 2,420 drone wings from the apiary in the settlement of Kuzemyn were classified into four clusters (Table 1). The classification correctness was 95.1%.
Mean values of indices for four clusters of drone wings (n=2420) from the apiary in the village of Kuzemin (Means SE; SE – standard error)
| Indexes | Cluster J1 (Means SE) | Cluster J2 (Means SE) | Cluster J3 (Means SE) | Cluster J4 (Means SE) |
|---|---|---|---|---|
| Ci | 2.121 0.016 | 1.473 0.010 | 1.754 0.017 | 1.572 0.012 |
| Dbi | 0.918 0.003 | 0.831 0.002 | 0.978 0.003 | 0.802 0.003 |
| Pci | 2.924 0.007 | 2.998 0.007 | 2.884 0.007 | 3.055 0.007 |
| Disc. sh. | 0.892 0.089 | 0.341 0.071 | 2.041 0.091 | -2.152 0.073 |
| Ri | 1.454 0.003 | 1.482 0.003 | 1.545 0.004 | 1.340 0.003 |
| Ci.3 | 1.827 0.008 | 1.678 0.005 | 1.900 0.009 | 1.823 0.007 |
| Ci.2.1 | 1.523 0.004 | 1.691 0.004 | 1.748 0.005 | 1.530 0.004 |
| C.2.2 | 4.179 0.035 | 4.449 0.027 | 3.795 0.030 | 4.182 0.027 |
| Ci.1.1 | 3.652 0.009 | 3.193 0.007 | 3.280 0.009 | 3.450 0.008 |
| Number of wings | 521 | 724 | 528 | 647 |
Coefficients of classification functions for the classification model obtained by classifying 2420 drone wings from the apiary in the village of Kuzemin into four clusters: J1, J2, J3, J4
| Indexes | Cluster J1 (p=0.21529) | Cluster J2 (p=0.29917) | Cluster J3 (p=0.21818) | Cluster J4 (p=0.26736) |
|---|---|---|---|---|
| Ci | 22.133 | 17.027 | 18.963 | 17.999 |
| Dbi | -0.644 | -4.039 | 19.372 | -7.571 |
| Pci | 58.796 | 61.022 | 58.030 | 63.565 |
| Disc. sh. | -6.210 | -6.232 | -5.993 | -6.765 |
| Ri | 258.367 | 254.261 | 260.643 | 244.915 |
| Ci.3 | 62.590 | 58.242 | 65.914 | 62.200 |
| Ci.2.1 | 164.664 | 173.363 | 175.484 | 163.616 |
| C.2.2 | -2.794 | -1.561 | -2.401 | -1.890 |
| Ci.1.1 | 123.166 | 113.561 | 113.813 | 118.017 |
| Constant | -697.280 | -664.160 | -704.640 | -662.470 |
The subspecies and population affiliation of wing phenotypes for the four clusters (J1, J2, J3, and J4) was determined using morphometric templates derived from previous studies conducted in 2022–2024 (Table 3).
Mahalanobis distances (MD) between the centroids of four clusters: J1, J2, J3, J4 and four drone wing templates obtained in previous studies
| Templates and clusters* | "UkrStep" | "HybrCaucas" | "Troiseck" | "Peschetz" | J1 | J2 | J3 | J4 |
|---|---|---|---|---|---|---|---|---|
| "UkrStep" | 0 | 3.52 | 2.76 | 1.84 | 2.90 | 1.23 | 2.66 | 3.21 |
| "HybrCaucas" | – | 0 | 3.18 | 3.92 | 3.57 | 2.72 | 4.61 | 1.45 |
| "Troiseck" | – | – | 0 | 2.41 | 1.38 | 2.74 | 2.57 | 2.78 |
| "Peschetz" | – | – | – | 0 | 2.86 | 1.83 | 1.71 | 3.56 |
| J1 | – | – | – | – | 0 | 2.98 | 2.61 | 2.71 |
| J2 | – | – | – | – | – | 0 | 2.73 | 2.44 |
| J3 | – | – | – | – | – | – | 0 | 3.87 |
| J4 | – | – | – | – | – | – | – | 0 |
Note: * – templates marked in transliteration: "Troiseck", "Peschetz" were obtained according to the work (Babenko et al., 2024b); "HybrCaucas", "UkrStep", – according to the work (Galatiuk et al., 2024b).
Four datasets, each comprising nine traits (indices) for clusters J1, J2, J3, and J4, were subsequently used as local templates (reference standards) for this apiary, designated in transliteration as "Troiseck(O)", "UkrStep(O)", "Peschetz(O)", and "HybrCaucas(O)", respectively. The symbol "O" is a specifying designation indicating that the experimental data originate from the eastern (Ost) part of Ukraine.
Drone wings for each studied colony are further classified into two clusters. The identification of each cluster was performed by comparing cluster centroids with the four local templates (Table 4).
Results of identification of 2,694 drone wings for 23 colonies from the apiary in the village of Kuzemyn. N denotes the number of wings in clusters, MD represents Mahalanobis distances between the centroids of clusters and templates, and DS represents the fractions of wings in clusters assigned to the templates: "Troiseck(O)", "Peschetz(O)", "HybrCaucas(O)", and "UkrStep(O)", respectively.
| Colony* | N wings | MD | Identification** | DS*** "Troiseck(O)", % | DS*** "UkrStep(O)", % | DS*** "Peschetz(О)", % | DS*** "HybrCaucas(О)", % |
|---|---|---|---|---|---|---|---|
| 305(1) | 68 | 1.59 | "Troiseck(O)" | 50 | 24 | 13 | 13 |
| 305(2) | 62 | 1.85 | "HybrCaucas(О)" | 23 | 0 | 0 | 77 |
| 307(1) | 57 | 2.24 | "HybrCaucas(О)" | 14 | 0 | 0 | 86 |
| 307(2) | 56 | 1.53 | "HybrCaucas(О)" | 7 | 21 | 0 | 71 |
| 310(1) | 53 | 2.06 | "Peschetz(О)" | 19 | 32 | 36 | 13 |
| 310(2) | 53 | 2.08 | "Peschetz(О)" | 4 | 9 | 87 | 0 |
| 313(1) | 62 | 1.28 | "HybrCaucas(О)" | 8 | 42 | 0 | 50 |
| 313(2) | 66 | 1.18 | "UkrStep(O)" | 0 | 94 | 6 | 0 |
| 314(1) | 55 | 1.61 | "HybrCaucas(О)" | 9 | 20 | 16 | 55 |
| 314(2) | 68 | >2.6 | Not identified | – | – | – | – |
| 316(1) | 52 | 1.73 | "UkrStep(O)" | 8 | 52 | 13 | 27 |
| 316(2) | 41 | >2.6 | Not identified | – | – | – | – |
| 322(1) | 71 | 1.78 | "Troiseck(O)" | 38 | 30 | 1 | 31 |
| 322(2) | 61 | 1.11 | "Troiseck(O)" | 80 | 5 | 15 | 0 |
| 324(1) | 64 | 1.65 | "HybrCaucas(О)" | 30 | 25 | 8 | 38 |
| 324(2) | 59 | 1.38 | "Troiseck(O)" | 90 | 0 | 5 | 5 |
| 329(1) | 60 | 2.28 | "Troiseck(O)" | 67 | 2 | 32 | 0 |
| 329(2) | 66 | 2.14 | "Troiseck(O)" | 45 | 14 | 8 | 33 |
| 331(1) | 67 | 1.31 | "HybrCaucas(О)" | 4 | 6 | 0 | 90 |
| 331(2) | 55 | 2.01 | "UkrStep(O)" | 38 | 36 | 20 | 5 |
| 332(1) | 56 | >2.6 | Not identified | – | – | – | – |
| 332(2) | 70 | 2.06 | "UkrStep(O)" | 1 | 54 | 37 | 7 |
| 339(1) | 68 | 0.99 | "UkrStep(O)" | 7 | 78 | 1 | 13 |
| 339(2) | 50 | 2.03 | "UkrStep(O)" | 0 | 68 | 32 | 0 |
| 312(1) | 52 | 1.28 | "Peschetz(О)" | 12 | 4 | 79 | 6 |
| 312(2) | 76 | 1.43 | "Peschetz(О)" | 0 | 14 | 78 | 8 |
| 342(1) | 69 | 2 | "Peschetz(О)" | 23 | 19 | 39 | 19 |
| 342(2) | 44 | 2.15 | "Troiseck(O)" | 77 | 2 | 20 | 0 |
| 317(1) | 23 | 1.95 | "HybrCaucas(О)" | 0 | 26 | 0 | 74 |
| 317(2) | 18 | 2.04 | "Peschetz(О)" | 6 | 61 | 17 | 17 |
| 330(1) | 61 | 2.26 | "Peschetz(О)" | 13 | 13 | 51 | 23 |
| 330(2) | 55 | 1.81 | "Peschetz(О)" | 2 | 9 | 89 | 0 |
| 338(1) | 55 | 1.52 | "HybrCaucas(О)" | 7 | 4 | 0 | 89 |
| 338(2) | 62 | 2.06 | "Troiseck(O)" | 40 | 45 | 5 | 10 |
| 326(1) | 92 | 1.78 | "HybrCaucas(О)" | 23 | 41 | 3 | 33 |
| 326(2) | 33 | 1.85 | "Troiseck(O)" | 55 | 21 | 24 | 0 |
| 375(1) | 59 | >2.6 | Not identified | – | – | – | – |
| 375(2) | 59 | 1.78 | "HybrCaucas(О)" | 5 | 51 | 0 | 44 |
| 327(1) | 53 | 1.87 | "HybrCaucas(О)" | 2 | 60 | 0 | 38 |
| 327(2) | 74 | 1.15 | "UkrStep(O)" | 3 | 88 | 9 | 0 |
| 315(1) | 54 | 1.61 | "UkrStep(O)" | 20 | 41 | 9 | 30 |
| 315(2) | 50 | >2.6 | Not identified | – | – | – | – |
| 309(1) | 70 | 1.91 | "UkrStep(O)" | 16 | 60 | 10 | 14 |
| 309(2) | 71 | 2.34 | "Peschetz(О)" | 4 | 15 | 80 | 0 |
| 333(1) | 60 | 2.01 | "UkrStep(O)" | 0 | 82 | 2 | 17 |
| 333(2) | 64 | 2.05 | "HybrCaucas(О)" | 27 | 16 | 11 | 47 |
Note: * – colony name, symbols 1 and 2 in parentheses indicate the cluster number; ** – name of the template to which this cluster is assigned; *** – wing classification was performed using a classifier obtained for the corrected array of 2420 drone wings.
The data presented in the “Identification” column indicate the assignment of cluster wing phenotypes to the “main breed,” while the “DS” columns show the proportions of wings that determine their conditional affiliation with the main and “introgressive” breeds (Table 4).
Presents summary statistics of the results of classification and identification of 46 clusters of drone wings. represents the Degree of Similarity of wing phenotypes in the clusters relative to the main "breed", expressed as the proportions of wings obtained using the four-cluster classifier
| Subspecific and population affiliation of clusters* | N**, (%) | N***, (%) () | N****, (%) () |
|---|---|---|---|
| A. m. carnica (Troiseck) | 9 (19.6) | 6 (13.04) | 3 (6.52) |
| A. m. carnica (Peschetz) | 8 (17.4) | 6 (13.04) | 2 (4.35) |
| Ukrainian steppe bees (UkrStep) | 11 (23.9) | 9 (19.57) | 2 (4.35) |
| A. m. caucasica (Hybrid caucasica) | 13 (28.3) | 8 (17.39) | 5 (10.87) |
| Not identified | 5 (10.9) | – | – |
| Together | 46 (100) | 29 (63.04) | 12 (26.02) |
Note: * – the assignment of wing cluster phenotypes to the main "breed" is determined according to Table 4; ** – the number of clusters assigned to the main subspecies or population; *** – the number of clusters where the share of wings of the main subspecies or population is ; **** – the number of clusters where the share of wings of the main subspecies or population takes values within: .
Presnets summary statistics of the assessment of hybridization of the main subspecies and populations by "admixtures" for the apiary in the village of Kuzemyn. represent the number of clusters for which the shares of "admixture" populations fall within the range of 10–50%
| Main population* | 1**, N1 | 2**, N2 | 3**, N3 |
|---|---|---|---|
| "Troiseck(O)" | "Peschetz(О)" 5 | "UkrStep(O)" 5 | "HybrCaucas(О)" 4 |
| "Peschetz(О)" | "Troiseck(O)" 4 | "UkrStep(O)" 5 | "HybrCaucas(О)" 3 |
| "UkrStep(O)" | "Troiseck(O)" 3 | "Peschetz(О)" 6 | "HybrCaucas(О)" 6 |
| "HybrCaucas(О)" | "Troiseck(O)" 5 | "Peschetz(О)" 2 | "UkrStep(O)" 13 |
Note: * – the assignment of wing cluster phenotypes to the main "breed" is determined according to Table 4; ** – 1, 2, 3 represent "admixture" populations.
Presents the results of classification of the combined array of drones wings (n = 2,420 — data from the present study; n = 74 — data for the template "Cc.307" from a previous study) into five clusters. The classification correctness was 94.4%
| Templates | Cluster 1 N, (%) | Cluster 2 N, (%) | Cluster 3 N, (%) | Cluster 4 N, (%) | Cluster 5 N, (%) | Number of wings, N |
|---|---|---|---|---|---|---|
| "Troiseck(О)" | 3 (0.58) | 1 (0.19) | 449 (86.18) | 9 (1.73) | 59 (11.32) | 521 |
| "UkrStep(О)" | 577 (79.70) | 0 | 0 | 72 (9.94) | 75 (10.36) | 724 |
| "Peschetz(О)" | 0 | 0 | 6 (1.14) | 278 (52.65) | 244 (46.21) | 528 |
| "HybrCaucas(О)" | 56 (8.66) | 480 (74.19) | 9 (1.39) | 0 | 102 (15.77) | 647 |
| "Cc.307"* | 13 (17.57) | 59 (79.73) | 1 (1.35) | 1 (1.35) | 0 | 74 |
| Together | 649 | 540 | 465 | 360 | 480 | 2494 |
Note: * – experimental data from the work (Galatiuk et al., 2023) were used.
Mahalanobis distances (MD) between the centroids of four patterns obtained in the current study and one in the previous study
| Templates | "Troiseck(О)" | "Peschetz(О)" | "UkrStep(О)" | "HybrCaucas(О)" | "Cc.307"* |
|---|---|---|---|---|---|
| "Troiseck(О)" | 0 | 2.86 | 3.16 | 2.75 | 3.53 |
| "Peschetz(О)" | – | 0 | 2.59 | 3.45 | 4.44 |
| "UkrStep(О)" | – | – | 0 | 2.48 | 2.81 |
| "HybrCaucas(О)" | – | – | – | 0 | 1.69 |
| "Cc.307"* | – | – | – | – | 0 |
Note: * – experimental data from the work (Galatiuk et al., 2023) were used.
For the identification of 22 clusters of drone wings from 11 colonies of an apiary located in the Dnipropetrovsk region, a classifier derived from the analysis of 2,420 wings from colonies in the settlement of Kuzemyn was applied (Table 2). The results of classification, phenotype identification for individual wing clusters, and the degree of probable hybridization are presented in Table 9.
Presents identification Results of identification of 988 drone wings from 11 colonies from an apiary in the Dnipropetrovsk region. N denotes the number of wings in clusters, MD represents Mahalanobis distances between the centroids of clusters and templates, and DS represents the fractions of wings in clusters assigned to the templates "Troiseck(O)", "UkrStep(O)", "Peschetz(O)", and "HybrCaucas(O)", respectively.
| Colony* | Wins, N | MD | Identific** | DS*** "Troiseck (O)", % | DS*** "UkrStep (O)", % | DS*** "Peschetz (О)", % | DS*** "HybrCaucas (О)", % |
|---|---|---|---|---|---|---|---|
| 14(IO)(1) | 24 | >2.60 | Not identified | – | – | – | – |
| 14(IO)(2) | 44 | 1.88 | "UkrStep(О)" | 2 | 98 | 0 | 0 |
| 32(IO)(1) | 33 | 1.83 | "HybrCaucas(О)" | 0 | 39 | 0 | 61 |
| 32(IO)(2) | 23 | 1.60 | "UkrStep(О)" | 13 | 74 | 13 | 0 |
| 75(IO)(1) | 40 | >2.60 | Not identified | – | – | – | – |
| 75(IO)(2) | 48 | 2.28 | "UkrStep(О)" | 15 | 58 | 19 | 8 |
| 44-32(1) | 34 | 2.52 | "UkrStep(О)" | 12 | 50 | 18 | 21 |
| 44-32(2) | 29 | 2.43 | "HybrCaucas(О)" | 0 | 3 | 0 | 97 |
| 66-55(1) | 46 | >2.60 | Not identified | – | – | – | – |
| 66-55(2) | 23 | 2.26 | "HybrCaucas(О)" | 0 | 22 | 0 | 78 |
| 67-32(1) | 56 | 2.04 | "UkrStep(О)" | 0 | 82 | 18 | 0 |
| 67-32(2) | 31 | 2.31 | "HybrCaucas(О)" | 0 | 58 | 0 | 42 |
| 69-55(1) | 31 | 2 | "Carnica"**** | – | – | – | – |
| 69-55(2) | 61 | >2.6 | Not identified | – | – | – | – |
| 46-32(1) | 42 | 1.62 | "HybrCaucas(О)" | 7 | 7 | 0 | 86 |
| 46-32(2) | 28 | 2.16 | "HybrCaucas(О)" | 0 | 4 | 0 | 96 |
| 11-32(1) | 54 | 1.83 | "UkrStep(О)" | 0 | 98 | 0 | 2 |
| 11-32(2) | 46 | >2.6 | Not identified | – | – | – | – |
| 27-14(1) | 60 | 2.01 | "HybrCaucas(О)" | 8 | 7 | 2 | 83 |
| 27-14(2) | 58 | >2.6 | Not identified | – | – | – | – |
| 31-32(1) | 92 | 1.83 | "HybrCaucas(О)" | 0 | 16 | 0 | 84 |
| 31-32(2) | 85 | 1.37 | "UkrStep(О)" | 12 | 69 | 4 | 15 |
Note: * – colony name, (IO) means that the queen was inseminated instrumentally, the symbols 1 and 2 in parentheses indicate the cluster number; ** – name of the template to which this cluster is assigned; *** – wing classification was performed using a classifier obtained for the corrected array of 2420 drone wings from the apiary in the village of Kuzemyne; **** – template was obtained according to data (Galatiuk et al., 2024b).
Summary statistics of the results of the assessment of hybridization of the main subspecies and populations by "admixtures" for an apiary located in the Dnipropetrovsk region
| Main population | 1*, N1 | 2*, N2 | 3*, N3 |
|---|---|---|---|
| UkrStep(O) | Troiseck(O) 1 | Peschetz(O) 3 | HybrCaucas(O) 2 |
| HybrCaucas(O) | Troiseck(O) 0 | Peschetz(O) 0 | UkrStep (O) 4 |
Note: * - 1, 2, 3 represent "admixture" populations; denote the number of clusters of "admixture" populations for which the shares of wings of "admixture" populations take values within the range of .

Shows the dendrogram of centroid distribution for the linear template "Cc.307"; local templates "Troiseck", "Peschetz", "UkrStep(Dn)"; and regional templates "Troiseck(О)", "UkrStep(О)", "Peschetz(О)", and "HybrCaucas(О)" based on Mahalanobis distances (full linkage rule). Points I and II represent critical points characterizing the distribution of templates within each group. The symbol (O) (Ost) indicates templates obtained for an apiary located in the eastern part of Ukraine, whereas (Dn) indicates templates obtained in the Dnipropetrovsk region. Templates "Cc.307" were taken from Galatiuk et al. (2023), and "Troiseck" and "Peschetz" from Babenko et al. (2024b).
Discussion
The study was based on experimental data obtained from a large apiary located in Kuzemyn, a border area of the Sumy and Poltava regions (Ukraine). The apiary has been continuously operating since 1991 and currently comprises approximately 860 colonies, of which about 450 are maintained exclusively for breeding purposes. Seasonal migration is performed within a radius of 3–6 km. The breeding program is focused on the maintenance, conservation, and reproduction of Ukrainian steppe honey bee populations. Annually, 5,000–9,000 queens are produced from 3–4 founder lines. The stationary queen-rearing unit includes approximately 220 colonies, and the nucleus system consists of about 2,450 mating units. Experimental groups are regularly established and maintained at a separate selection apiary (up to 200 colonies), where phenotypic traits are recorded over a two-year period. In the third year, selected groups are used to form the breeding core for the subsequent propagation cycle. Mating of queens for selection purposes is conducted under controlled conditions on isolated islands in the Dnipro River near the settlements of Zhovnyne and Kuchuhury. The surrounding area (within a radius of 5–10 km) contains approximately 700 additional colonies from other apiaries.
A total of 2,964 drone wings from 23 colonies were analyzed to assess the population structure of the apiary. At the preliminary stage, wings were classified into four clusters using discriminant analysis. The resulting sets of morphometric indices were considered provisional local templates (reference standards) for the studied apiary. In parallel, wings from each colony were classified into two clusters. Cluster identification was performed by comparing cluster centroids with the provisional templates using Mahalanobis distances (MD). Clusters with MD≥2.6 relative to all templates were excluded from further analysis. After filtering, 2,420 wings were retained. At the main stage, the filtered dataset (n=2,420) was subjected to repeated classification into three- and four-cluster solutions. Classification correctness was 96.5% and 95.1%, respectively. Cluster identification was performed using morphometric templates obtained in previous studies (Babenko et al., 2024b; Galatiuk et al., 2023, 2024) (Tables 1 and 3). In the three-cluster solution, two clusters were assigned to Apis mellifera carnica, while the third cluster corresponded to the Ukrainian bee population (UkrBee). In the four-cluster solution (J1–J4), the structure of the two A. m. carnica clusters remained largely unchanged. However, the cluster previously assigned to UkrBee split into two distinct components, allowing a more refined identification of Ukrainian bee populations. The clusters were interpreted as follows: J1 — *A. m. carnica *(Troiseck line), J3 — A. m. carnica (Peschetz line), J2 — Ukrainian steppe population (UkrStep), and J4 — hybrid population (HybrCaucas) (Tables 1 and 3).
Taking into account the conclusions of previous studies (Galatiuk et al., 2024b) regarding the possible influence of the subspecies A. m. caucasica on wing phenotype, this assumption was tested for populations identified by the "HybrCaucas(O)" and "UkrStep(O)" templates (Table 8). Additionally, to evaluate potential introgression, a combined dataset including 2,420 wings and 74 wings from the reference template "Cc.307" (originating from a queen from Georgia, presumably of the Megrelian population) was subjected to extended classification analysis. The optimal solution corresponded to a five-cluster model (Table 7). It was established that 74.2% of the wings of the "HybrCaucas(O)" template and 79.7% of the "Cc.307" wings belonged to the same cluster, indicating high morphological similarity. This is confirmed by the low Mahalanobis distance between their centroids (MD=1.69). At the same time, the distance between the centroids of "Cc.307" and "UkrBee" is 1.60 (Galatiuk et al., 2024b), which provides grounds to assume the presence of genetic introgression of the subspecies A. m. caucasica into the UkrBee and HybrCaucas(O) populations, particularly in the genomic component associated with wing morphology. Accordingly, cluster J4 can reasonably be denoted as HybrCaucas (Caucasian hybrid). Additionally, cluster 5 should be noted, in which substantial proportions of wings of the templates "Peschetz(O)" (46.21%), "HybrCaucas(O)" (15.77%), "Troiseck(O)" (11.32%), and "UkrStep(O)" (10.36%) were identified. This indicates partial admixture of phenotypes of the subspecies *A. m. carnica *with Ukrainian populations.
As noted previously (Babenko et al., 2024), further classification of wings for individual colonies within the already identified two clusters into a larger number of clusters lacks clear biological justification. However, its use is appropriate as an analytical tool for assessing the structure of phenotypic variability. Using a classifier derived from the full dataset, each cluster was re-evaluated in terms of proportional assignment to the four local templates. These results (Tables 4 and 9) provide quantitative estimates of phenotypic composition and admixture. Given the high density of surrounding colonies and the historical context of the apiary, the co-occurrence of A. m. carnica, A. m. caucasica, and populations of the Ukrainian steppe bee is non-random and systematic. Consequently, ongoing hybridization among these groups is expected and is consistently detectable at the phenotypic level using morphometric analysis of drone wings.
Analysis of the data presented in Table 4 indicates that, for 29 of the 41 studied drone wing clusters, the majority of wings (>50%) were assigned to a primary template, which determined the overall cluster identification. A smaller proportion of wings was assigned to secondary templates. For seven clusters, two wing groups were observed, with more than 20% of wings assigned to each of two different templates, one of which corresponded to the primary template. In five clusters, wings were approximately evenly distributed among three templates. It is important to note that all drone wings within a given cluster originate from the same queen, i.e., they are produced by the same type of haploid genome and are therefore phenotypically similar. Consequently, the observed distribution across multiple templates should be interpreted as an estimate of genomic hybridization. The proportion of wings assigned to non-primary templates (DS) is interpreted as the degree of hybridization (DH) with other subspecies or populations, following the approach of Babenko et al. (2024a). For 11 clusters, the primary template proportion (DS>80%) indicated no detectable hybridization. Eleven clusters exhibited hybridization with one other subspecies or population. Thirteen clusters showed hybridization with two other populations. Five clusters involved three contributing populations, with minor components having DS>10% (Table 4). Exceptions were observed for clusters no. 327(1) and no. 67-32(2) (Table 9), where Mahalanobis distances (MD=1.87 and 2.31) assigned wings to the "HybrCaucas(O)" and "UkrStep(O)" templates, respectively, while the majority of wings were classified under the opposite template according to the discriminant model. Additional classification for 46 clusters showed that, in 26% of cases, the proportion of wings assigned to the primary population ranged between 30% and 50%, which is insufficient for confident interpretation (Table 5). These findings highlight the complex interactions among populations and the limitations of the method for precise determination of hybridization type in some cases.
The proportion of clusters assigned to A. m. carnica versus Ukrainian populations (UkrStep(O) and HybrCaucas(O)) favored the latter, with a ratio of 17:24 (41.5% vs. 58.5%; Table 4). Clusters associated with "Troiseck(O)" and "Peschetz(O)" exhibited comparable levels of hybridization. In contrast, hybridization among UkrStep(O) and HybrCaucas(O) populations predominated over hybridization involving A. m. carnica, suggesting that introgression from
*A. m. carnica *has not yet exceeded 50%.
In the Dnipropetrovsk apiary, only one cluster was assigned to the "Carnica" population (Carpathian ecotype; Table 9). Seven clusters (31.8%) were identified as UkrStep(O), and another seven (31.8%) as HybrCaucas(O). Among UkrStep(O) clusters, six showed 15%–25% contributions from introgression components, while four HybrCaucas(O) clusters exhibited similar introgression contributions; in one case, the introgression component exceeded 25%. Overall, interactions among Ukrainian bee populations predominated (6 vs. 4 cases involving introgression components of A. m. carnica; Table 10). Due to limited sample sizes and the fact that 27.3% of clusters remained unidentified, these results are insufficient for a quantitative assessment of population proportions within the small geographic range of the Dnipropetrovsk apiary. Nevertheless, Ukrainian populations clearly predominate over A. m. carnica in this region.
The availability of specific morphometric templates for individual lineages of the subspecies A. m. carnica and A. m. caucasica enables analysis of drone wing morphology within a two-level hierarchical framework: population level → subspecies level. This is clearly illustrated by a dendrogram constructed for two apiaries, reflecting relationships both among populations and among subspecies. The dendrogram was built using morphometric templates obtained in both the present and previous studies. The templates were divided into two main groups: (1) A. m. carnica: "Troiseck(O)", "Troiseck", "Peschetz", "Peschetz(O)";
(2) A. m. caucasica/Ukrainian bees: "Cc.307"+"HybrCaucas(O)" corresponding to* A. m. caucasica*, and "UkrStep(O)"+"UkrStep(Dn)" presumably corresponding to Ukrainian bees (A. m. macedonica). This grouping is consistent with the discriminant analysis performed through additional classification of the combined dataset (Table 7). However, it creates difficulties in determining the subspecies affiliation of the templates "UkrStep(O)" and "UkrStep(Dn)", which likely belong to populations of Ukrainian steppe bees, subspecies
*A. m. macedonica *(Galatiuk et al., 2024). The Mahalanobis distance between the subgroups ["Cc.307"+"HybrCaucas(O)"] and ["UkrStep(O)"+"UkrStep(Dn)"] is approximately 3.8, indicating a substantial difference. The divergence between groups (1) and (2) is even greater (MD≈4.8), which explains the higher phenotypic similarity of Ukrainian steppe bees to A. m. caucasica (evolutionary lineage O) rather than to A. m. carnica (evolutionary lineage C). Thus, at present, there are insufficient grounds for reliable determination of the subspecies affiliation of the populations designated as UkrStep(O) and UkrStep(Dn).
To identify bees of the subspecies A. m. caucasica, the authors used experimental data from a single queen (no. Cc.307) as a reference template, which may raise concerns regarding the relevance of this approach and the reliability of the results obtained. However, the validity of this approach is supported by data reported in Janashia et al. (2025), indicating that: “Georgia may offer the most viable mainland strategy for conserving pure A. m. caucasica populations, as large-scale introductions of non-native subspecies have been minimal. Hypothetically, the exclusive use of native honeybees by Georgian beekeepers supports the expectation that the country still harbours A. m.caucasica populations in their purest form”. The aforementioned work also notes that: “… the natural distribution of A. m. caucasica remains poorly defined due to insufficient data from key regions, including Georgia, Türkiye, Armenia, Azerbaijan, and Russia (Northern Caucasus)”. In light of the theory of *Apis mellifera *dispersal (Han et al., 2012), Ukraine may also be considered within this broader dispersal range. According to morphometric and genetic data, 29 subspecies of Apis mellifera are divided into four evolutionary lineages: A, M, C, and O. Lineages C (A. m. carnica) and O (A. m. caucasica) share a common origin, forming a distinct root branch within the species; however, the exact stage of divergence between them has not yet been established. Despite the existence of different scenarios for the eastward expansion of these lineages into Europe, there are grounds to assume possible shared dispersal routes, which would facilitate hybridization between the subspecies A. m. carnica and
A. m. caucasica.
The results of this and previous studies (Yarovets et al., 2023; Galatiuk et al., 2024a, b) indicate that apiaries in the Lviv, Zhytomyr, Kyiv, Kharkiv, Poltava, and Dnipropetrovsk regions of Ukraine host multiple Apis mellifera populations attributed to the subspecies A. m. carnica, A. m. caucasica, as well as local Ukrainian populations (UkrStep, UkrBee) and their hybrids (HybrCaucas). No apiaries with isolated purebred populations were identified, indicating the absence of a clearly defined range of purebred populations for these “breeds.” The factors determining the current breed composition of bees in the studied apiaries are complex. The geographic landscape from the Caucasus to the Carpathians does not present significant barriers to honey bee dispersal in northern, northwestern, and western directions, including along the coasts of the Azov and Black Seas, which is consistent with dispersal theory (Han et al., 2012). Additionally, during the 1960s–1980s, hundreds of thousands of
A. m. caucasica queens were introduced into Ukraine as part of state-supported agricultural programs, which substantially complicated breeding efforts and contributed to the widespread formation of hybrid populations. Given these circumstances, several strategies for breeding program implementation can be recommended (Uzunov et al., 2022a; Brascamp et al., 2024). Where isolated mating sites are available, the establishment of closed “purebred” populations is advisable. In open populations, breeding efforts should focus on hybrids that have a quantitative advantage over other “breeds.” In both scenarios, instrumental insemination should be applied, as it is an effective method for stabilizing apiary structure and improving economically important traits within a relatively short period (Uzunov et al., 2022b; Du et al., 2023, 2024).
Conclusion
The analysis of a large sample of drone wings from the apiary in Kuzemyne enabled the development of robust morphometric templates for four drone wing phenotypes, which can be regarded as local reference standards. The application of these templates for identifying wing morphology in the geographically distant Dnipropetrovsk region demonstrated their suitability for use across a broader spatial context.
The authors propose the use of these templates as morphometric markers alongside genetic markers in breeding programs. The analysis of the obtained results, together with the proposed approach for assessing the degree of hybridization of the main breed through introgressed components, provides evidence for a substantial presence of the subspecies A. m. caucasica in the eastern and central regions of Ukraine.
Adherence to ethical principles. The requirements of the recommendations of the EU Directive 2010/63/EU on the protection of animals used for scientific purposes have been met.