1 Introduction

The honey bee Apis mellifera, a world-wide distributed species, is an insect with enormous environmental and socio-economic importance. Various hypotheses place the origin of this species in Africa or the Middle East (Ruttner 1978; Whitfield et al. 2006; Weinstock et al. 2006; Cridland et al. 2017) with a subsequent diversification into lineages, subspecies, and ecotypes to cope with environmental challenges throughout most ecosystems of the world. Ruttner (1988) identified four major lineages: African (A), North Mediterranean (C), Middle Eastern (O), and West European (M). Using molecular methods, Franck et al. (2001) later discerned a fifth lineage (Y) based on Ethiopian honey bees. However, there are controversies in grouping of subspecies to an evolutionary lineage using molecular methods (De la Rúa et al. 2009). According to Arias and Sheppard (1996), subspecies from Northeastern Africa are more related to Middle Eastern groups of subspecies, forming the lineage O with a partial disagreement with preceding classical morphometric classification (Ruttner 1988). In Northeastern Africa and the Middle East, there are several contact zones between lineages A, O, and Y (Cridland et al. 2017), because the lineage A represents African bees (Ruttner 1988; Whitfield et al. 2006), Y exists in Ethiopia, and the Middle Eastern lineage O extends up to Egypt, Sudan, and Somalia (Franck et al. 2001; El-Niweiri and Moritz 2008). However, there was no indication of lineage Y in Sudan (El-Niweiri and Moritz 2008), which shares an extended border with Ethiopia. Later, Alburaki et al. (2011) classified A. m. syriaca as Z subgroup of African lineage A along with A. m. lamarckii and Alburaki et al. (2013) designated A. m. syriaca as a separate lineage Z differing from lineages A and O using microsatellite analysis. More recently, Cridland et al. (2017) showed five global lineages and confirmed that lineage Y exists in close proximity with A and O lineages by analyzing A. m. jementica samples from Yemen and Saudi Arabia to represent the Ethiopian lineage Y.

Ethiopia owns an estimated ten million honey bee colonies (Girma 1998), of which about six million are managed (Central Statistical Agency [CSA] 2017) by two million smallholder farmers (Anand and Sisay 2011). Considering the distribution and classification of honey bee subspecies in Ethiopia, several reports exist that are partly inconsistent (Radloff and Hepburn 1997; Amssalu et al. 2004; Meixner et al. 2011). Amssalu et al. (2004) reported five subspecies including A. m. scutellata, A. m. jementica, A. m. monticola, A. m. bandasii, and A. m. weyi-Gambella. Radloff and Hepburn (1997) described three morphologically distinct clusters that correspond to A. m. jementica, A. m. bandasii, and A. m. scutellata. However, Radloff and Hepburn (2000) separated Ethiopian honey bees from scutellata, and Meixner et al. (2011) described them as a distinct subspecies A. m. simensis. These contrasting reports were merely based on classical morphometric analysis.

In order to elucidate the evolutionary history and level of differentiation among honey bees, a variety of different approaches have been used. The methods range from standard morphometric and behavioral analyses (Ruttner 1988), mitochondrial (e.g., Garnery et al. 1993), and microsatellite DNA comparisons (e.g., Franck et al. 2001) to wing geometric morphometrics (Nawrocka et al. 2018) and genome-wide SNP analyses (Whitfield et al. 2006; Zayed and Whitfield 2008; Wallberg et al. 2017). Forewing venation, which can be efficiently measured (Tofilski 2004; Nawrocka et al. 2018), carry sufficient information from both parental lines to discriminate different groups of honey bees (Meixner et al. 2013). On the other hand, the intergenic region between cytochrome oxidase subunit I and subunit II (COI-COII) in the mitochondrial DNA of honey bees was proven to be informative to elucidate evolutionary history of maternal lineages, simple, and widely documented in various population studies.

From an evolutionary point of view, Ethiopia, which is located in Northeast Africa at close proximity to the Middle East, is an important region for honey bees. The country is known for its diverse agroecological zones (AEZs), ranging from humid to arid, and elevations from > 100 m below sea level to 4500 m above sea level (masl). In view of the general consensus on the presence of multiple lineages and subspecies in the surroundings of the Ethiopian plateau, two hypotheses can be discussed. First, the country may be inhabited by multiple subspecies as reported by Amssalu et al. (2004) that belong to both lineages A and O. The alternative hypothesis would be that the Ethiopian honey bees may not constitute more than one subspecies, but rather form a unique subspecies A. m. simensis, as proposed by Meixner et al. (2011), that belongs to the evolutionary lineage Y (Franck et al. 2001).

Based on the situations listed, there is a need for a comprehensive study on the Ethiopian honey bees in order to establish a solid basis on their evolutionary relatedness, diversity, and geographic distribution; paving ways for genetic conservation and apicultural development.

In the present study, we explored the diversity of Ethiopian honey bees using wing geometric morphometrics and COI-COII mitochondrial DNA analyses. Wing venation is less influenced by environmental factors and can efficiently separate honey bee populations (Meixner et al. 2013), whereas the intergenic region of COI-COII varies among honey bee lineages in sequence length and frequency of characteristic P and Q motifs (Garnery et al. 1993; Meixner et al. 2013). Therefore, evolutionary lineages, subspecies, and mitochondrial haplotypes were investigated and their distribution was assessed, based on AEZs and geographic location associated with the samples. The in-depth analysis of honey bee phylogeny in this region may further contribute to the design of sustainable conservation and breeding strategies. These may become necessary since human activities, such as habitat modification (Härtel and Steffan-Dewenter 2014), selective breeding, introduction, and replacement of native populations (Michener 1975; Meixner et al. 2010), as well as colony marketing between agroecological zones (Teweldemedhn and Yayneshet 2014) impact honey bee populations.

2 Material and methods

2.1 Site selection and sampling

Samples were collected in two main geographic regions of Ethiopia, Tigray regional state (north) and Wendogenet local area (south), representing three and two AEZs, respectively (Figure 1). Tigray is a regional state that is located in the northern part of the Federal democratic republic of Ethiopia whereas Wendogenet local area is a small geographic area located in the southern part of the country sharing areas between the regional states of Oromia and South Nations Nationalities and Peoples (SNNP). Within the Tigray, three local areas (Mugulat, Werie, Koyetsa), each extending from highland (≥ 2100 masl) and midland (1800–2100 masl) to lowland (< 1800 masl) AEZs, were selected, whereas sampling in Wendogenet local area was limited to midland and lowland AEZs (Figure 1) due to political instability and insecurity during sample collection in 2018. Therefore, nine sampling sites were selected out of highland (Mugulat, Kolageble and Koyetsa), midland (Adikebero, Tsedya, Adidaero), and lowland (Werie river nearby the main bridge, Simret, Aditsetser) AEZs from three local areas in Tigray regional state and two sampling sites out of midland (Wendogenet) and lowland (Entaye) AEZs from Wendogenet local area (Table I; Figure 1). Local areas are named after one of the sample sites it contained that better represent the area (Table I). A distance of at least 10 km was maintained between sampling sites to minimize interference of colony dispersal through swarming and mating. By selecting six colonies from each site and 10 worker bees per colony, a total of 660 samples were collected out of 66 colonies (hereafter denoted as this study) from the eleven sites. Acquiring export permit from Ethiopian Institute of Biodiversity in accordance with national proclamation 482/2006 (Federal parliament 2006), samples were imported to the University of Hohenheim (Germany) for genetic and morphometric analyses.

Figure 1.
figure 1

Study site maps: A total of nine sample sites were selected to represent three agroecological zones (AEZs)-highland, midland, and lowland-in Tigray region from northern Ethiopia (upper), indicated by red, yellow, and green markers, respectively. Two sample sites were included from the southern part, represented by Wendogenet (midland) and Entaye (lowland).

Full size image
Table I Description of the sample colonies used in this study (number sampled and codes given), sample site location (administrative zone, latitude, longitude), elevation, agroecological zone (AEZ), temperature, and precipitation
Full size table

Tigray was selected based on its beekeeping potential (Central statistical agency [CSA] 2017), indication of two subspecies (Amssalu et al. 2004), and its geographic importance (proximity to habitat of the lineage O) while Wendogenet was included based on its distance from the north, its mountain system, and proximity to areas of the lineage A. Tigray in particular and Ethiopia in general composes of mainly three AEZs (highland, midland, lowland). Average annual temperature ranges 15.28 to 18.9°C in the highlands, 17.4 to 20.5°C in the midlands, and 20.8 to 25.39°C in the lowlands of Tigray. Average annual precipitation is 467 to 508 mm in the highlands, 461 to 495 mm in the midlands, and 547 to 620 mm in the lowlands of Tigray (Global agroecological zones). Moreover, the highland study areas are mainly mountains covered with dry montane vegetation including Juniperus procera, Olia africana, and Becium grandiflorum; whereas the midlands are covered mainly with Acacia itbica, Otostegia integrifolia, Cordia africana, and Acacia albida. The lowlands have small hills mainly covered with deciduous thorny shrubs and trees of acacia species, evergreen large riparian fig tree species, and Ziziphus abyssinica (Figure S1).

2.2 Wing geometric morphometric analysis

The right forewings of the samples (n = 660) were detached and digitalized for geometric morphometric analyses while the thorax of each bee was labeled and stored in 80% ethanol at – 20 °C for molecular analysis. Images of wings were taken using ZEISS Axiocam 105, which was mounted on a stereomicroscope, and processed with ZEN lite 2.1 software (Zeiss, Germany). Nineteen landmarks were manually measured by one person on each image using Identifly 1.1.0 (Tofilski 2017). Each measurement was double-checked for consistency. Landmark coordinates were superimposed using full procrustes fit that involves translation, rotation, and scaling of landmarks using MorphoJ (Klingenberg 2011). This allows to conduct multivariate analysis on shape (Zelditch et al. 2004), besides to univariate analysis on centroid size. Wings shape was described by configuration of nineteen landmarks, and wing size by centroid size. Both aligned coordinates and centroid size were averaged within colonies and the later statistical analysis was based on the average values. Samples representing four major lineages (16 M, 37 C, 86 A, and 50 O) and eleven African subspecies were included from Nawrocka et al. (2018) as references in lineage and subspecies classification, respectively. In addition, 140 images of lineage Y representing fourteen colonies of A. m. simensis (Meixner et al. 2011) were obtained from Oberursel Bee Research Institute, Germany.

2.3 Mitochondrial DNA analysis

Total DNA was extracted from the thorax of 66 worker honey bees representing 66 colonies using a modified CTAB protocol (Rusterholz et al. 2015). Primer sequences E2 and H2 were compared with established sequences of different subspecies available at the National Center for Biotechnology Information (NCBI) in order to amplify the COI-COII intergenic region by polymerase chain reactions (PCR). Because of consistent mismatches with the reference sequences, both forward and reverse primers were modified (E2i:5′-GGCAGAATAAGTGCATTG-3′, H2i:5′-CAATATCATTGATGACCA-3′). PCR reactions were performed at the following conditions: initial denaturation at 95 °C for 120 s; followed by 35 cycles of denaturation at 95 °C for 30 s, 48 °C annealing for 30 s, and a two-step elongation at 60 °C for 15 s and 65 °C for 90 s; and final elongation at 65 °C for 5 min. Afterwards, an aliquot of the PCR amplicons were digested with the restriction endonuclease DraI (37 °C overnight) and separated by agarose gel (4%) electrophoresis. PCR amplicons were purified by ethanol precipitation and sent for sequencing (Eurofins Genomics). The forward primer and a longer variant of the reverse primer (H2i-seq: 5′-CAATATCATTGATGACCAATTG-3′) were used for sequencing. Restriction fragment bands (Figure S2) and amplicon lengths were estimated with the gel images and cross-checked by in silico digesting the sequences with DraI (Table II).

Table II Description of haplotypes based on DraI restriction, PQ combination, total size (bp: base pair) of COI-COII intergenic region, and % (out of 62 included samples)
Full size table

Using the progressive alignment algorithm in CLC Main Workbench 7.6.4 (QIAGEN, Aarhus Denmark), a total of 83 sequences were aligned on nucleotide level, of which 62 were generated from the samples of this study (four samples failed amplifications) and 21 were reference sequences representing evolutionary lineages A, C, M, O, and Y from NCBI (Figure 2, Table S1). Phylogenetic analysis was conducted by maximum likelihood, for which we have initially identified the nucleotide substitution model that best fit to the data using the Model Test option implemented in MEGA X (Kumar et al. 2018). Based on the lowest Bayesian Information Criterion (BIC) score, Tamura-Nei model (Tamura and Nei 1993) was used to infer evolutionary history. A bootstrap consensus tree was inferred from 500 replicates to represent the evolutionary history (Felsenstein 1985) by applying neighbor-joining and BioNJ algorithms using all sites. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates were collapsed. Phylogenetic analysis was carried out using MEGA X (Kumar et al. 2018).

Figure 2.
figure 2

Maximum likelihood evolutionary analysis conducted on 83 nucleotide sequences of COI-COII using Tamura-Nei model (Tamura and Nei 1993). Numbers at the nodes represent bootstrap values. Analysis was conducted in MEGA X (Kumar et al. 2018). Sample bees of this study are represented by four-digit identification numbers. The first digit (right) stands for individual bee sample (one of 10 bees), the second digit represents colony numbers 1 to 6 in every site, the third digit stands for the AEZ of the specific site (1 refers to highland, 2 to midland, 3 to lowland), and the fourth digit stands for local area: 1 for Mugulat, 2 for Werie, 3 for Koyetsa, and 4 for Wendogenet. The reference samples are labeled by accession number and haplotype code/subspecies name/country of origin as found in the repository.

Full size image

2.4 Statistical analysis

Sample colonies were classified into lineages and subspecies using canonical variate analysis (CVA) based on wings shape. The relationship of centroid size with elevation (masl), latitude (decimal degree), and longitude (decimal degree) was evaluated using Pearson correlations. Afterwards, linear regression analyses were used to verify the dependence of centroid size on the factors correlated with it. The data were tested for normal distribution using goodness of fit. In addition, multivariate analysis of variance (MANOVA) was run to compare this study samples with reference lineages and subspecies based on wings shape.

Statistical significance of variation in the distribution of subspecies and haplotypes was tested using nominal logistic regression model considering elevation (masl), latitude (decimal degree), longitude (decimal degree), and local area (four local areas) as factors. Furthermore, variation in the subspecies and haplotype distribution among AEZs (three AEZs) and sampling sites (eleven sites) was assessed using the same model. The distribution of haplotypes along AEZs and local areas were plotted in a graph including a 95% ellipse. All statistical analyses were performed using Statistica 13® (from Tibco Data Science) and JMP® Pro 15 (from SAS).

3 Results

3.1 Geometric morphometric analysis

Since it was suggested that Ethiopian honey bees constitute a separate lineage Y (Franck et al. 2001), firstly, we tested if lineage affiliation could be confirmed by geometric morphometric analysis. Thus, the reference samples for A. m. simensis were compared with reference samples of lineages A, C, M, and O. Canonical variate analysis based on forewing shape clearly differentiated A. m. simensis from the other four lineages (Figure 3). The analysis allowed to correctly classify, with cross-validation, all reference colonies of A. m. simensis. Moreover, Mahalanobis distance between A. m. simensis and lineage A (46.28) was greater than the Mahalanobis distance between lineages A and O (38.20) (Table IV). Hence, A. m. simensis is further referred to as lineage Y. Forewing venation differed significantly between the five evolutionary lineages (MANOVA: Wilks Λ  0.0006, F = 26.2, P < 0.001). In pairwise comparisons, all lineages differed markedly from each other (Table IV). Secondly, the samples from this study were projected into the space of the canonical variates of the lineages that were obtained earlier. In this case, our samples overlapped mostly with lineage Y (Figure 3). Thus, it can be concluded that the sample was most similar to lineage Y followed by lineage A, based on Mahalanobis distances (Table IV). Specifically, 91% of the sample colonies were classified as lineage Y, while the rest was classified as lineage A (Table S2). In the next step, the samples from this study were compared with the reference samples of African honey bee subspecies. A .m. simensis was most similar (88%) to the samples from this study (Table S3, Table III; Figure 4); however, some colonies were classified as A. m. scutellata (9%), A. m. monticola (1.5%), and A. m. unicolor (1.5%). None of the factors among elevation, latitude, longitude, local area, sampling sites, and AEZs showed a significant influence (P > 0.05) on the distribution of subspecies. Hence, samples classified as A. m. simensis were abundantly present in all AEZs and local areas, whereas those classified as non-simensis subspecies were sporadic (Table S3).

Figure 3.
figure 3

Discrimination of honey bee lineages using first and second (a) or second and third (b) canonical variate (CV). A 95% confidence level ellipse plot of honey bee lineages showing the separation of this study samples from evolutionary lineages A, C, M, and O but overlapped with lineage Y. The discrimination of lineage Y is particularly well visible on the graph of the second and third canonical variates. Each marker represents one colony.

Full size image
Table III Squared Mahalanobis distances between this study sample and references of eleven African honey bee subspecies (lower triangle) and statistical significance of pairwise comparisons between the groups (upper triangle)
Full size table
Figure 4.
figure 4

Canonical variate analysis on forewing shape: A 95% confidence level ellipse plot of honey bee subspecies discriminated this study samples from reference subspecies of A. m. scutellata, A. m. monticola, A. m. adansonii, and A. m. jemenitica but mostly (88%) overlapped with A. m. simensis. Each marker represents one colony.

Full size image

Forewing size positively correlated with increasing elevation (r = 0.64; P < 0.01) and longitude (r = 0.39; P < 0.01) but it was not significantly correlated with latitude (r = 0.04; P = 0.77). The centroid size depended on elevation as confirmed by a linear regression analysis (F = 43.7; P < 0.01; R2adj. = 0.4); however, longitude could explain little of this association (R2adj = 0.14).

3.2 Analysis of COI-COII intergenic region

Our genetic analyses on the COI-COII intergenic region provided new insights into haplotype diversity, distribution, and evolutionary lineage of Ethiopian honey bees.

We identified five mitochondrial haplotypes that included four previously reported, Y (Y1 and Y2), A (A1), and O (O5’) and one new haplotype that we have denoted as Y3. Haplotype Y2 accounted for majority (62.9%), Y3 was the second most abundant (19.3%) of the sample whereas O5’ (1.6%) stood the least. There was a significant difference in the distribution of these haplotypes between the local areas (X2 = 32.1, P < 0.01) but not between elevations, geographic coordinates, sampling sites, and AEZs (P > 0.05). All samples that belonged to the haplotypes Y1, Y2, and Y3 lay within a 95% ellipse, whereas the rest fell outside. Haplotype Y2 was found throughout the sampling sites, but A1 was observed in two highland AEZs (Kolageble and Koyetsa) and Y1 was located mainly along a river valley (Werie) in the central zone of Tigray region as well as the lowland of Wendogenet area. On the other hand, haplotype Y3 extended throughout all elevations in Koyetsa area of Northwest Tigray. Although most sampling sites had one or two haplotypes, the highland of Koyetsa had four (A1, Y1, Y2, Y3) and the midland of Werie had three (Y1, Y2, Y3) haplotypes (Figure S3). Based on single nucleotide polymorphism (SNP), the five haplotypes were resolved into 33 variants (Table II).

Phylogenetic analysis using maximum likelihood method showed that the majority (74.2%) of the Ethiopian samples clustered with lineage Y. However, some (17.7%) colonies mainly from Koyetsa local area were grouped with lineage O. Furthermore, 8.1% grouped with haplotypes of lineage A (Figure 2). Looking at the P sequence of COI-COII, 79.3% of the samples had this segment with a size ranging from 60 to 66 bp. This sequence occurred zero to three times in the COI-COII region among the samples. Four (6.4%) samples had this part of the sequence with 67 bp (P0) but 14.3% (9) lacked it (Table II). Overall, the P and Q combinations in this study were PQ (58.7%), PQQ (20.6%), P0Q (4.8%), P0QQ (1.6%), Q (3.2%), and QQ (11.1%); considering 67 bp as P0 and 60–66 bp as P. Depending on the frequency and combination of the P and Q sequences, COI-COII PCR amplicons exhibited three band lengths ranging from 637 to 1025 bp (data not shown). Samples with the supposedly ancestral (P0) sequence were distributed among haplotypes Y1, Y2, and Y3.

4 Discussion

Based on geometric morphometrics and mitochondrial DNA analysis, our comprehensive study provides novel insights into the position of Ethiopian honey bees within the evolutionary lineages (Figure 3, Figure 2) and neighboring subspecies (Figure 4). Four evolutionary lineages of A. mellifera were postulated by Ruttner (1988), based on classical morphometric analysis that categorized all African honey bee samples of that time as lineage A. Later, Franck et al. (2001), who conducted analysis of COI-COII mitochondrial DNA intergenic region, reported that Ethiopian honey bees represent a new linage (Y) based on consistent deletion of some nucleotides at a specific position in the P sequence (termed as P2) differing from west European lineage M (P) and sub-lineage A of Atlantic coast (P1). The colonies sampled during the present study predominantly belong to evolutionary lineage Y, and thereby support the findings of Franck et al. (2001). Moreover, we noticed that A. m. simensis markedly differs morphologically from other subspecies of lineage A (Table III; Figure 4). Lineage Y shows a slightly increased level of differentiation to lineage A compared to lineage O (Table IV) and separates clearly morphologically (Figure 3a).

Table IV Squared Mahalanobis distances between this study samples and references of honey bee lineages A, C, M, O, and Y (lower triangle) and statistical significance of pairwise comparisons between the groups (upper triangle)
Full size table

Our data on wing morphometrics confirmed an earlier study (Meixner et al. 2011) that honey bees of Ethiopia belong to a unique subspecies denoted as A. m. simensis whereas some of our samples were classified as other subspecies (A. m. scutellata, A. m. monticola). Earlier reports indicated up to five subspecies in the country, including A. m. scutellata, A. m. jementica, and A. m. monticola (Amssalu et al. 2004). Similarly, Radloff and Hepburn (1997) described three clusters: A. m. jementica (North), A. m. bandasii (Central), and A. m. scutellata (South). The different results could be related to methodological variations. For example, Amssalu et al. (2004) considered the presence of clusters as a separate subspecies and reported five subspecies whereas Meixner et al. (2011) compared distances between clusters with reference subspecies although both studies have similar coverage and distribution of samples in the country. A. m. scutellata is distributed in Savannah areas of Eastern and Southern Africa, whereas A. m. monticola is found in distributed patches on East African mountains (Smith 1961; Ruttner 1988; Radloff and Hepburn 2000; Gruber et al. 2013). Therefore, it could be predicted that A. m. simensis would be similar to one or more of these subspecies. Interestingly, similarity of A. m. simensis to these subspecies is relatively low (Table III, Table S3) and examined specimens of this study formed a single cluster (Figure 3, Figure 4). This suggests that our knowledge about the distribution of honey bee subspecies in this geographic region is still incomplete. It is important to mention that the current version of IdeniFly software has small size of reference samples (Tofilski 2017).

Considerable proportion of the samples lacked the P sequence within the analyzed mtDNA-fragment. Complete deletion of the P sequence was reported to be a characteristic of the Southeastern European honey bee lineage C, while a full length P0 sequence (67 bp) is typical for most of the African and Middle Eastern honey bees (Garnery et al. 1993).The absence of a P sequence was also reported from honey bees in neighboring Sudan (El-Niweiri and Moritz 2008), West African Benin (Amakpe et al. 2018), and haplotype Y2 of Ethiopia (Franck et al. 2001). A small fraction of the samples from this study was found to have a P0 sequence, without a clear regional or ecological distribution pattern. In contrast to morphometrics, the mitochondrial study indicated that some of our samples from the Koyetsa area are closer to haplotypes of lineage O, which extends from the Middle East to the neighboring Northeast African countries (El-Niweiri and Moritz 2008; Alattal et al. 2014).

Dissimilarity between Ethiopian honey bees and neighboring populations can be explained by geographic and climatic isolation of the Ethiopian plateau, unique agroecological features within the country, and the pastoral farming system along the boundaries in the South, East, and West. Ethiopia is generally classified as central-northern highlands and peripheral arid lowland rangelands (FAO 1986). The lowlands are characterized by pastoral livestock production system (Coppock 1994), where beekeeping is less favored. Strips of suitable habitats along river valleys are probably too narrow for substantial gene flow, so the presence of lineage A among the Ethiopian samples of this study was negligible. In contrast, the national border in the north does not impair the gene flow of bees, and is characterized by relatively advanced beekeeping practices such as use of frame hives and colony splitting and marketing (Nuru 2002; Teweldemedhn and Yayneshet 2014). Hence, a considerable proportion (56%) of the sample in northwest zone of Tigray clustered with lineage O in a phylogenetic tree analysis (Figure 2), indicating an overlap of the lineage O and Y habitats.

These findings strengthen the hypothesis that Ethiopia is situated on one of the potential spreading routes of honey bee lineages between Africa and the Middle East. However, Mahalanobis distance based on forewing venation indicated that lineages Y and A are closer to each other than to O (Table IV), suggesting that the origin and distribution routes of the honey bee need further investigation. Cridland et al. (2017) supported the origin of A. mellifera to be in Northeast Africa or the Middle East, with the ancestral population giving rise to A and Y lineages. However, the precise placement of the two lineages is not clear. The origin of the samples that they have used (Harpur et al. 2014) differs from those used in describing the lineage Y (Franck et al. 2001). A recently conducted review (Dogantzis and Zayed 2019) indicated lineage Y to be of Asian origin in contrast to the initial description of the lineage Y to represent Ethiopian honey bees (Franck et al. 2001). The tendency to associate lineage Y with the Middle East could be arisen due to the fact that Ethiopian honey bees were assumed to be A. m. jementica when characterizing the lineage Y, although later it was renamed as A. m. simensis after being discriminated from A. m. jementica and other neighboring subspecies (Meixner et al. 2011). Hence, a further investigation will have to clarify if lineage Y extends beyond Ethiopia, be it Middle East or neighboring African countries. So far, there was no indication of lineage Y in Sudan (El-Niweiri and Moritz 2008), which shares extended border with Ethiopia. Some confusion arose among different studies, resulting into ambiguities when classifying subspecies such as A. m. jementica, A. m. lamarckii, and A. m. syriaca into evolutionary lineages. Ruttner (1988) showed that A. m. jementica, which was classified into African lineage A along with A. m. lamarckii, occupies several countries in Africa (Sudan, Chad, Somalia) and Middle East (Saudi Arabia, Yemen, Oman). The lineage O was reported to be extended up to Egypt, Sudan, and Somalia (Franck et al. 2001; El-Niweiri and Moritz 2008) among African countries; Yemen and Saudi Arabia (Alattal et al. 2014) as well as into ranges of A. m. syriaca (Ruttner 1988; Wallberg et al. 2014) in the Middle East. Contrastingly, Alburaki et al. (2011) classified A. m. syriaca (honey bee samples from Syria, Lebanon, and Iraq) together with A. m. lamarckii as Z subgroup of African lineage A; and later, Alburaki et al. (2013) designated A. m. syriaca as a separate lineage Z.

Our analysis of forewing venation across AEZs provided evidence of local differentiation, indicative as signs for local adaptation of these honey bees. The results showed that samples from the highland have substantially larger forewings than those from lowland areas, which is consistent with previous findings on body size variation at different elevations in Ethiopia (Meixner et al. 2011). Local adaptation to high elevations has been shown for east African mountain bees, denoted as A. m. monticola, which showed behavioral and morphological differences compared to the adjacent lowland populations (Gruber et al. 2013). Interestingly, remarkable genetic differentiation, mainly restricted to only 1.4% of the genome, has been described for these bees (Wallberg et al. 2017). Therefore, this ecological inclination seen in Ethiopian bees will be object of further studies to elucidate potential nuclear genetic differentiation among them.

Our study revealed that Ethiopia harbors several haplotypes, which seem to be influenced by human activities in terms of their distribution. This was evident in local areas that are the main destinations of traded colonies (Werie and Koyetsa) being populated by several haplotypes in contrast to local areas in the eastern zone that mainly supply colonies (Figure S3). However, this study focused in Tigray region of northern Ethiopia with a purposeful inclusion of two sampling sites from the southern part of the country. Honey bee colony marketing involving different AEZs is a common practice in the northern part of the country (Nuru 2002; Teweldemedhn and Yayneshet 2014). Considering the unique mating nature that enables genetic exchange over wide distances (Jensen et al. 2005) and the highly mobile behavior of African honey bees (McNally and Schneider 1992), it is also likely that there could be extensive gene flow and hybridization. This is possible in light of the local patchy forests that cannot support isolated populations as standalone habitats; allowing gene flow between AEZs as discussed by Gruber et al. (2013). The connection of habitats has been enhanced by decades-long extensive work on ecosystem restoration in Ethiopia (Balehegn et al. 2019). Generally, African honey bees are characterized by low genetic differentiation, which could be resulted from their migratory behavior (Franck et al. 2001; Fuller et al. 2015), besides anthropogenic influences.

In conclusion, the results of this study supported the hypothesis that Ethiopian honey bees are distinct both at lineage and subspecies levels, despite significant morphometric variability and diverse mitochondrial haplotypes. Further research on nuclear DNA will provide deeper insights into the level of hybridization and potential local adaptation in the Ethiopian honey bees.