Genetic diversity and novel haplotypes of Apis mellifera jemenitica on the Arabian Peninsula: insights from mtDNA markers

Abstract

The genetic diversity of Apis mellifera jemenitica populations collected from the Arabian Peninsula (Saudi Arabia, Yemen, and Oman), Jordan, and Ethiopia, was examined using three mtDNA markers: 1- Cytochrome b (Cyt b), 2- Cytochrome c oxidase I (COI) and 3- The intergenic region located between the cytochrome c oxidase I & II (COI-COII). DNA was extracted from 44 samples, amplified for each region using classic PCR, and the resulting amplicons were sequenced using Sanger technology at both ends. Sequences were verified and aligned, and Maximum-Likelihood phylogenetic analyses were conducted with reference sequences from other subspecies. The in silico DraI mtDNA COI-COII (DmCC) test was applied to the COI-COII sequences to identify evolutionary lineages and haplotypes. Moreover, COI-COII haplotype network analyses were conducted to assess the intra- and inter-genetic relationships between samples and references. Based on the Cyt b marker, most samples cluster within the African lineage (A) near lamarckii and syriaca (Sub-lineage Z) subspecies. Few samples from Ethiopia and Yemen were closely related to simensis and scutellata clades. The COI gene separated jemenitica samples (Bootstrap = 97) from subspecies of other lineages (C and O). The DmCC test revealed a P0Q2 structure in the intergenic region for all samples, with a distinct 18 bp deletion in the P0 element observed in two Ethiopian and one Yemeni samples, suggesting litorea or simensis origin. A total of 13 COI-COII haplotypes were identified, among which 8 haplotypes were novel: Saudi Arabia (1), Yemen (3), Oman (1), and Ethiopia (3), with a haplotype diversity (H) of 0.980. Furthermore, molecular-variance parsimony in COI-COII confirmed a distant genetic relationship between Ethiopian samples versus samples of the Arabian Peninsula. The haplotype network analysis suggests a higher intra-jemenitica diversity than previously understood with a syriaca ancestry to this clade. These findings offer crucial insights into the conservation of A. m. jemenitica and its role in preserving biodiversity in arid ecosystems. Additionally, the data enhance our understanding of the genetic diversity of A. m. jemenitica and its evolutionary connections with other neighboring African subspecies.

1 Introduction

Apis mellifera jemenitica Ruttner (1976) is a major honey bee subspecies of the African evolutionary lineage (A). In his book “Biogeography and Taxonomy of Honeybees,” Ruttner discussed this poorly studied subspecies under the Chapter “Honeybees of Tropical Africa” alongside a handful of other African subspecies such as lamarckii, litorea, scutellata, adansonii, monticola, capensis and unicolor (Ruttner, 1988). Based on morphometrical characteristics, jemenitica is considered one of the most intriguing and significant subspecies of A. mellifera, notable for its extreme morphology and ecology. It was not recognized as a taxonomic unit until 1976 (Ruttner, 1988). This subspecies extends into a large, hot, arid territory of 4,500 km from East to West, covering many countries, including Yemen, Saudi Arabia, Oman, Ethiopia, Chad, Somalia, and Sudan (Dutton et al., 1981; Ruttner, 1988; Ruttner et al., 1978). Its distribution is strongly associated with the availability of flowering plants and appropriate nesting sites. This subspecies has been used in the Arabian Peninsula for beekeeping since 2,000 BA (Alqarni et al., 2011). While still poorly studied, jemenitica is described as specifically adapted to the ecological conditions of Yemen and surrounding regions (Al-Ghamdi and Nuru, 2013b). Phenotypically, jemenitica is among the smallest honey bee subspecies, typically displaying black or dark brown coloration, which may represent an adaptation for conserving energy and resources in arid environments (Ruttner, 1988). Other studies reported jemenitica as almost the size of Apis cerana (Oldroyd and Wongsiri, 2006) with a yellow abdomen and grey to brown bands (Al-Ghamdi and Nuru, 2013b). This suggests potential wide phenotypical and genetic diversity within jemenitica populations. Behaviorally, jemenitica is known for its relatively calm nature, making it easier to manage during beekeeping activities (Al-Ghamdi and Nuru, 2013b). Like African subspecies, jemenitica exhibits a tendency to swarm, which is a natural reproductive strategy. Jemenitica bees are efficient foragers, capable of sourcing nectar and pollen from various plants, including those that thrive in arid conditions. Their foraging behavior is adapted to maximize resource utilization in environments with limited floral availability (Al-Ghamdi and Nuru, 2013a). This subspecies has developed several physiological adaptations to cope with high temperatures, such as more efficient thermoregulation and behaviors that minimize exposure to heat. The foraging quality and efficiency of jemenitica in its native range of distribution make it a vital pollinator for many native plants, contributing to local biodiversity and the health of ecosystems.

Beekeeping practices in the regions of jemenitica often incorporate traditional methods passed down through generations. It is estimated, at least for the jemenitica populations of Saudi Arabia, that 70% or more are kept in traditional hives, often made of palm tree trunks, clay pots, and simple mud (Al-Ghamdi and Nuru, 2013b). Like many honey bee species, jemenitica faces threats from habitat loss, pesticide use, climate change, diseases, and the importation of foreign subspecies to its native repartition areas. Protecting natural habitats, promoting sustainable agricultural practices, and raising awareness of the importance of this subspecies are crucial for its conservation.

Honey bee diversity worldwide was first evaluated using morphological traits, including tergite pigmentation, body and wing dimensions, pilosity, tomentum width, cubital index, and proboscis length (Ruttner, 1988; Ftayeh et al., 1994; Meixner et al., 2007; Ruttner et al., 1978; Kauhausen-Keller et al., 1997; Kandemir et al., 2000; Arias et al., 2006; Alburaki and Alburaki, 2008). Four evolutionary lineages of honey bees were identified based on morphological traits: the West Mediterranean lineage (M), North Mediterranean lineage (C), African lineage (A), and Oriental lineage (O) (Ruttner, 1988; Ruttner et al., 1978; Rinderer, 1986). Although some recent studies have reported discrepancies, all four morphometrical lineages have been molecularly confirmed through mitochondrial and microsatellite markers (Smith and Brown, 1988; Smith and Brown, 1990; Cornuet et al., 1991; Garnery et al., 1992; Garnery et al., 1993; Estoup et al., 1995; Moritz et al., 1994; Garnery et al., 1998a; Garnery et al., 1998b; Meixner et al., 2000; Smith et al., 1991; Bouga et al., 2005). There are still unresolved discrepancies between mitochondrial DNA (mtDNA) and nuclear DNA analyses, particularly concerning the O and Y lineages (Ilyasov et al., 2020). For instance, the Y lineage has been reported to extend into the southern Arabian Peninsula and the Horn of Africa (Dogantzis and Zayed, 2019). However, this lineage has not been reported using morphometric traits or mtDNA markers (Ruttner, 1988; Ruttner et al., 1978; Kandemir et al., 2006). Currently, 33 honey bee subspecies, including jemenitica, have been identified in their native regions across Europe, Western Asia, and Africa (Sheppard and Meixner, 2003; Ruttner, 1988; Chen et al., 2016; Ilyasov et al., 2020).

The mitochondrial DNA (mtDNA) is a valuable tool in phylogeography and genetic studies due to its low or no recombination, relatively high evolutionary rate, conserved structure, and uniparental heredity (Avise et al., 1987; Harrison, 1989; Avise, 2000). The cytochrome b (Cyt b) gene is an essential component of the mitochondrial genome in bees and other organisms (Toparslan et al., 2020; Tobe et al., 2010). This gene encodes a protein critical in the electron transport chain, vital for cellular respiration and energy production in the mitochondria. The cyt b gene is often studied for its genetic variability among honey bee lineages and subspecies (Crozier and Crozier, 1992; Alburaki et al., 2011; Pinto et al., 2004). Its variation can help us understand honey bee evolutionary relationships, population structure, and local adaptation to specific environments. Variations in the cyt b gene can also offer insights into the genetic health of honey bee populations by assessing the impact of environmental stressors and diseases on bee populations. Coupled with other markers, the cyt b gene can serve as a valuable tool in bee research, aiding in understanding the genetics, evolution, and conservation of honey bees. Similar to the cytochrome b (cyt b), cytochrome c oxidase I (COI) is a mitochondrial coding gene commonly used to investigate the genetic diversity of honey bee subspecies (Kandemir et al., 2006).

The DraI mtDNA COI-COII test, referred to as (DmCC), is a widely used mitochondrial test in honey bee genetics (Garnery et al., 1993; Smith and Brown, 1990; Cornuet et al., 1991). This test is conducted on the mtDNA intergenic region situated between both cytochrome c oxidase I and II. This region’s polymorphisms can determine, to a great extent, the evolutionary lineage, subspecies, and haplotype, rendering this test one of the most comprehensive and efficient tests to explore honey bee genetic diversity (Alburaki et al., 2011; Rortais et al., 2011; Arias and Sheppard, 1996; De la Rua et al., 1998; Franck et al., 2001; Palmer et al., 2000; Loucif-Ayad et al., 2014; Franck et al., 2000; Hall and Smith, 1991; Abrahamovich et al., 2007; El-Niweiri and Moritz, 2008).

The genetic diversity of jemenitica has marginally been explored during the last decades (Alghamdi et al., 2012; Alattal et al., 2014b). This study aims to evaluate the genetic diversity of jemenitica using mitochondrial DNA markers and to assess its evolutionary relationships with neighboring subspecies. We hypothesize that jemenitica populations exhibit significant genetic diversity due to their broad geographic range and unique adaptations to various arid environments. Our data identified novel haplotypes, indicated structural differences between jemenitica and other neighboring subspecies, and a wider jemenitica genetic diversity within the African lineage than previously presumed.

2 Materials and methods

2.1 Honey bee samples

A total of 44 honey bee samples were sampled from five different countries: 1- Saudi Arabia (n = 22), Yemen (10), Ethiopia (7), Oman (3), and Jordan (2), Figure 1; Table 1. The sample size of 44 individuals from five countries was limited by logistical constraints but deemed sufficient to capture a broad representation of jemenitica populations across the study regions. One honey bee worker was sampled per colony. The region and city of each sample are detailed in Table 1. All samples originated from small-scale traditional apiaries, mainly belonging to stationary beekeepers. All samples (n = 44) were analyzed for cyt b, and 15 samples were studied for the COI gene and the intergenic region COI-COII. Worker bee samples were sampled in absolute ethanol and stored at 20°C at the Bee Research Unit at King Saud University for further downstream analyses. Due to mid-experiment technical difficulties and a change of institutions by the principal investigator, only 15 samples of the original 44 samples were available to run for COI and COI-COII regions.

FIGURE 1

2.2 DNA extraction

Genomic DNA was extracted from the thoraxes of honey bee workers using the DNeasy^® Blood and Tissue Kit (Qiagen, Germany), following the manufacturer’s protocol. Briefly, thoraxes were dissected using a scalpel and forceps and kept in 2 mL microcentrifuge tubes. Two glass beads and 180 μL ATL buffer were added to the tubes, and tissues were homogenized for 5 min at speed 10 (Equivalent to 6,000 cycles/minute) in a Bullet Blender (Next Advance, Inc., Averill Park, NY). A 20 μL of Proteinase K (10 mg/mL) was added per sample, vortexed, and incubated at 56°C overnight. Tube contents were transferred to DNeasy mini spin column for a DNA column extraction according to the manufacturer’s protocol. DNA extracts were stored at −20°C for further use.

2.3 Amplification of the coding genes cyt b and COI

All samples (n = 44) were amplified at both ends for the cyt b gene using the universal honey bee set of primers: Cytb-F: 5′- TAT​GTA​CTA​CCA​TGA​GGA​CAA​ATA​TC-3′ and Cytb-R: 5′-ATT​ACA​CCT​CCT​AAT​TTA​TTA​GGA​AT-3′ (Crozier et al., 1991; de Oliveira Francisco et al., 2001; Sheppard et al., 1994; Alburaki et al., 2011). Fifteen samples among the 44 were amplified for the COI gene using the following set of primers: COI-F: 5′- TCT​ATA​CCA​CGA​CGT​TAT​TC-3′ and COI-R: 5′- GAT​CAA​TAT​CAT​TGA​TGA​CC -3′ (Smith et al., 1997). The PCR reaction was conducted in a 12 µL reaction composed of 6 µL BioRad Master Mix, 0.5 µL of each primer (10 µM), and 5 µL of nuclease-free water. The PCR cycling for both genes was 95°C for 3 min, followed by 35 cycles of (95°C for 10 s, 55°C for 30 s) and a permanent hold at 4°C. The standard two-step PCR amplification was efficient and produced optimal bands, eliminating the need for a final extension phase. Negative controls were included in all PCR reactions to ensure the absence of contamination, and positive controls were used to verify amplification success. PCR products were verified by electrophoresis on a 1.5% agarose gel.

2.4 Amplification of the intergenic non-coding region COI-COII

The intergenic non-coding region of the mtDNA located between the cytochrome c oxidase I and II was studied on the same fifteen samples evaluated for the COI gene. This region was amplified at both ends using the traditional set of primers: E2-F: 5′-GGC​AGA​ATA​AGT​GCA​TTG-3′, H2-R: 5′-CAA​TAT​CAT​TGA​TGA​CC-3′ (Alburaki et al., 2011; Rortais et al., 2011; Achou et al., 2015; Meixner et al., 2000; Garnery et al., 1993; Franck et al., 2001). The PCR reaction for this set of primers was carried out in a 26 µL reaction composed of 12.5 µL of Promega’s 1X Taq buffer, 1.9 mM MgCl2, 0.2 mM of each dNTP, 0.14 µM of each primer, 16.5 µL of nuclease-free water, 1.25 U/reaction of Taq, and 1.4 µL of DNA. PCR cycling consisted of the following parameters: 95°C for 2 min, followed by 30 cycles of (95°C for 30 s, 54°C for 30 s, 72°C for 30 s), and 72°C extension phase for 10 min. Negative controls were included in all PCR reactions to ensure the absence of contamination. Similar to both previous genes, amplicons were run on 1.5% agarose gel for size and quality checks.

2.5 Amplicon sequencing and alignment

PCR products were purified using Applied Biosystems ExoSAP-IT kit and sequenced at both ends using an automated sequencer (Applied Biosystems 3730XL DNA Analyzer, United States). Sequencing quality was assessed using (3730XL Data Collection Software 5), with ambiguous bases trimmed and only high-quality reads included in the analysis. Each sample’s forward and reverse sequences for all three regions (Cyt b, COI, and COI-COII) were paired after reversing the reverse sequence, aligned, and manually checked using Geneious Prime^® 2024.0.7 (Drummond et al., 2011). Nucleotide ambiguities at both sequence ends were trimmed, and all sequences generated from this study (n = 74) were deposited at the National Center for Biotechnology Information (NCBI) GenBank.

2.6 Data analysis

Geographical mappings of the populations’ distributions were conducted using the Tableau Public platform (https://www.tableau.com). Sequence alignment, verification, generation of FASTA files, and phylogenetic analysis were carried out using Geneious Prime^® 2024.0.7 (https://www.geneious.com). Sequences of both coding genes (Cyt b and COI) were separately aligned by MUSCLE standard alignment algorithm (Edgar, 2004), and reconfirmed by Geneious alignment using a Cost Matrix of 65% similarity (5.0/-4.0) with a “Global alignment with free end gaps” option, then trimmed at both ends to eliminate non-identified nucleotides. However, the COI-COII non-coding sequences (n = 15) were first sorted out according to the genetic structure and length of the intergenic region (Q, P, P0, PQ_x, P0Q_x), in which (x) is a tandem repeat of the Q fragment, and analyzed following the steps described in previous studies (Madella et al., 2020; Alburaki et al., 2023). All sequences were blasted against the NCBI databases, and the highest blasting scores were reported in this study, along with their reference sequences. Novel COI-COII haplotypes were named according to the universal nomenclature detailed in a recent study (Alburaki et al., 2023).

Generation of haplotype networks and computation of nucleotide diversity and Tajima’s statistic were conducted in the R environment (R Core Team, 2022) (v. 2024.09.0) using molecular-variance parsimony available in the libraries: “ape” and “pegas” (Paradis, 2010; Paradis and Schliep, 2019). Maximum-likelihood phylogenetic trees were conducted using the Tamura-Nei distance model (Tamura and Nei, 1993) employing the “PHYML” plugin (Guindon et al., 2010) with default parameters: Number of bootstraps = 100, automatic estimation of transition ratio, proportion of invariable sites, and Gamma distribution parameter, with a default number of substitution rate categories of four. This study prioritized the maximum likelihood method for constructing phylogenetic trees, given the limited number of sequences available and the assumption of varying evolutionary rates among the samples (Zou et al., 2024). When relevant, genetic trees were rooted using A. cerana sequences (FJ229480) or genetic clades from distant non-African lineages. Only NCBI-verified and accurate reference sequences of other subspecies were used in this study, including previously published jemenitica haplotypes from the same studied regions (Supplementary Table S1). For some genes, the absence of neighboring reference sequences (e.g., simensis, litorea, monticola) has restricted the scope of our analyses.

3 Results

3.1 Phylogenetic analysis of cyt b gene

All forty-four studied samples from five countries (Figure 1) were amplified for the cyt b gene with a 359 bp length after trimming ambiguous nucleotides at both ends, Table 1. Out of 44 sequences, 35 blasted at 100% with reference samples representing three different subspecies: jemenitica, simensis, and lamarckii. The nine other samples blasted at 99.7% as follows: 1- Jordan (n = 2): syriaca, 2- Ethiopia (n = 1): lamarckii, 3- Saudi Arabia (n = 1): iberiensis, and 4- Yemen (n = 5): jemenitica, simensis, unicolor, lamarckii, Table 1. A total of 11 cyt b sequence groups were recorded, characterizing six different subspecies, Table 1. In the phylogenetic analysis, most samples clustered in an independent African lineage (A) clade with lamarckii and separated from the African sub-lineage (Z) with a 36-bootstrap value, Figure 2A. Samples from Jordan, the native region of syriaca, clustered with syriaca references (Figure 2A). Two Yemeni samples (Yem4, 8) and three Ethiopian samples (Eht4, 5, 6) showed closer genetic proximity to simensis and scutellata, Figure 2A. The intra-sample phylogenetic tree, rooted with A. cerana separated (33-bootstrap) the Jordanian samples (Jor1, 2) from the rest as well as Omani samples and one single Ethiopian sample (Eth1) with a 64-bootstrap value, Figure 2B. The sample (Sau4), which clustered within the M lineage (Figure 2A), remained quite distant (90-bootstrap) from other samples, Figure 2B. All forty-four cyt b sequences are available at NCBI GenBank under accession numbers PQ460982 to PQ461025, Table 1.

FIGURE 2

3.2 COI sequence analysis

Fifteen samples were sequenced for the COI coding gene, Table 2. The verified and retained sequence length of the COI gene was 218 bp. NCBI-verified reference sequences in African subspecies for the COI gene were much less abundant, which restricted the scope of this analysis. Five sequence groups were recorded within the studied samples, Table 2. All sequences except one (Sau15) blasted at 100% with reference sequences of various subspecies: jemenitica, syriaca, lamarckii, simensis, scutellata, monticola, capensis, mellifera and ligustica, Table 2. Sau15 exhibited a SNP (C replacing T) at the gene end (155 bp). Based on the sequences of this gene, Yem4 and Eth5 are the most genetically distant among the studied samples, Figure 3A. Eth4, however, formed an independent cluster with a bootstrap value of 34. In contrast, other samples formed a solid clade (97-bootstrap) with a jemenitica reference (OQ348419) away from other subspecies belonging to the C, O, and M lineages, Figure 3A. In the intra-sample phylogenetic analysis, Eth4 and Eth5 were the most distant samples, followed by Yem4 (28-bootstrap), Yem2 (93-bootstrap), and Sau14 (75-bootstrap), while the rest of the samples exhibited high genetic similarity with weak bootstrap values, Figure 3B. All fifteen COI sequences are available at NCBI GenBank under accession numbers PQ461026 to PQ461040, Table 2.

FIGURE 3

3.3 COI-COII structure, haplotype and phylogeny

The same fifteen samples studied for COI were sequenced for their intergenic COI-COII regions. The COI-COII sequences were captured, varying in length from 818 to 839 bp, Table 3. All sequences exhibited one of the typical African COI-COII structures (P0Q2). Interestingly, a deletion of 18 pb was identified in the P0 element of three samples (Eth4, Eth5, Yem4), Table 3. This distinct and significant deletion matches those found in reference jemenitica sequences, all collected from Ethiopia (FJ477998, FJ477999, FJ478000, FJ478001, FJ478002, FJ478003), Table 3. A total of 13 COI-COII haplotypes were identified, among which eight were novel, Table 3. Six samples fully blasted (100%) with previously identified jemenitica haplotypes from Saudi Arabia (KSA4e, KSA6b, KSA6d, KSA4f) and one African haplotype from the United States (A1-818-6-USA), Table 3. The highest blasting scores (99.4%–99.9%) for the novel haplotypes were all with jemenitica references and one sample (Yem8) with syriaca, Table 3. Only reference sequences of the same COI-COII structure (P0Q2) were selected and run in the phylogenetic analysis. The studied samples clustered into two distinct clades supported by a 100-bootstrap value, Figure 4A. Group 1 comprised most samples clustering with reference sequences of jemenitica, syriaca, and Group 2, representing previously recorded Ethiopian haplotypes. The latter branch was separated from other North African subspecies and supported with a 30-bootstrap value, Figure 4A. The intra-sample COI-COII variation is higher than what was seen for the two other genes. Eth4 and Eth5 were the most distant haplotypes, followed by Yem4 (99-bootstrap) and Sau14 and Sau13 (Same haplotype; 100-bootstrap), Figure 4B. The novel Omani haplotype (A1-838-6-OMN) showed greater similarity to Eth1, another novel haplotype (A1-839-6-ETH) identified in this study (Figure 4B). All fifteen COI-COII sequences are available at NCBI GenBank under accession numbers PQ461041 to PQ461055, Table 3.

FIGURE 4

3.4 Haplotype diversity and network

A significant number of COI-COII haplotypes (n = 13) were identified in a relatively small population of 15 samples. The haplotype diversity (H) was 0.980, consistent with the corrected (H), as no non-African COI-COII haplotypes were identified (Table 3; Figure 5). The nucleotide diversity (π) within the studied samples was (0.025) with a Tajima’s statistic (D = −1.25, p > 0.05) indicating no deviation from neutrality or significant selection activity, Figure 5. Interestingly, the COI-COII haplotype network suggests a significant genetic distance between Ethiopian haplotypes (Eth1, Eth4, Eth5) and the rest, Figure 5. All other haplotypes seem to have diverged from each other through (≥1) SNPs, Figure 5. Within this cluster, Yem2 occupied a central position from which other haplotypes may have diverged, Figure 5.

FIGURE 5

An assorted number of (P0Q2) NCBI reference sequences (n = 47, Supplementary Table S1) from various subspecies were selected to study the inter-genetic diversity of our samples. The COI-COII haplotype network distinguished four groups of haplotypes: 1- Haplotypes of the M lineage (M4p, M4k, M4m) known to characterize both mellifera and iberiensis subspecies, 2- Eastern African haplotypes (Y1a, Y1b, Y1c, Y1d, Y2a) sampled from Ethiopia with which clustered three of the studied samples (Eth4, Eth5, Yem4), 3- North African haplotypes (A4p, A9i, A9l) characterizing intermissa, including two African haplotypes (A3-830-4-USA, A1-829-4-USA) recently identified in the United States, and 4- The majority of the Arabian Peninsula’s samples including those of Jordan clustering with reference samples of syriaca and previously identified jemenitica haplotypes from Saudi Arabia, Figure 6. This analysis identified a scutellata haplotype (A4) sequence (FJ477987) as a potential shared ancestry from which all other groups deviated, Figure 6.

FIGURE 6

4 Discussion

Based on 24 morphometrical characteristics, discriminant analysis sharply distinguished jemenitica samples of Saudi Arabia from other neighboring subspecies such as litorea, lamarckii, and syriaca (Alattal et al., 2014a). All samples studied from Saudi Arabia (n = 179) exhibited the P0 elements in their COI-COII intergenic region (Alattal et al., 2014a), which makes them genuine members of the African lineage (Alburaki et al., 2011). Similar findings were recorded for jemenitica samples (n = 16) from the isolated Island of Socotra (Alattal et al., 2019). In the light of our COI-COII results, all analyzed samples (n = 15) exhibited the P0 element and thus belong to the African lineage (A), blasting with jemenitica and syriaca reference sequences, Table 3. The structural differences in the sequence following the tRNA_Leu gene have facilitated lineage and haplotype differentiation in A. mellifera. Specifically, samples from the (C) lineage and its haplotypes lack the P fragment, whereas the (A) lineage includes the P0 fragment (68 bp), and the (M) lineage contains the 54 bp P fragment (Cornuet et al., 1991; Garnery et al., 1993; Franck et al., 2000; Garnery et al., 1998a; Garnery et al., 1995; Rortais et al., 2011; Alburaki, 2011). The intertwining between jemenitica and syriaca observed in this study is yet to be understood by further investigations. Recent mtDNA sequencing studies have pointed out this genetic relationship between jemenitica, syriaca, and lamarckii (Boardman et al., 2020; Alghamdi and Alattal, 2020). Phylogenetic analysis of the cyt b offers different nuances. Most samples clustered within the African clade of lamarckii separated with a weak bootstrap value of 36 from the syriaca sub-lineage “Z” (Alburaki et al., 2013). This result indicates an undoubted genetic relationship with the African subspecies syriaca but still closer proximity with the lamarckii subspecies of Egypt, Figure 2. Only three Ethiopian and two Yemeni samples deviated from this cluster, showing closer relationships with simensis and scutellata,Figure 2. Notably, three of these five samples (Eth4, Eth5, Yem4) carried a significant deletion of 18 bp in their COI-COII P0 element, Table 3. Previous investigations on three Apis m. jemenitica samples from Ethiopia identified similar deletions in the P0 element (Alattal et al., 2019; Franck et al., 2001). It is safe to conclude that these deviating samples, at least for the Ethiopian samples (Eth4, 5, 6), are not jemenitica but simensis. This conclusion is supported by a 100% blast of their cyt b sequences with simensis references, Table 1. In contrast, the three other Ethiopian samples (Eth1, 2, 3) represent jemenitica, Table 1; Figure 2. Yemeni samples (Yem4, 8) showed closer similarity with unicolor and simensis than jemenitica, Table 1; Figure 2. From a biogeographical perspective, such phenomena are often observed in zones of subspecies interactions. The three major neighboring subspecies of jemenitica across the Red Sea in East Africa are simensis, monticola, and litorea (Ruttner, 1988). Unfortunately, very few sequences of the two latter subspecies were available to shed more light on the overlapping or interaction of these African subspecies with jemenitica. Even when available, most sequences were partial and could not be retained for analysis.

The COI gene offers a similar outcome, with the three most distant samples (Yem4, Eth4, 5) from the jemenitica cluster, a cluster that was firmly (Bootstrap = 97) separated from subspecies of all other lineages (C, M, and O), Figure 3. Nonetheless, at sequence blasting level, the COI gene seems to provide weak discriminative power as sequences fully blasted with multiple subspecies, Table 2. However, this may have occurred due to the relatively short length (218 bp) of our sequences in the COI gene. The COI phylogeny seems to fully agree with the cyt b, clustering the same four samples (Eth4, 5, Yem2, 4) as the most distant sequences from the rest, even at the intra-sample level, Figure 4. Interestingly, this observation was recorded to a great extend (Eth4, 5, Yem4) for the COI-COII intergenic region, Figure 4B. The COI-COII analysis showed a clear separation of four major clades: 1- M lineage, the African lineage (A) subdivided to 2- jemenitica (Group 1) which groups the majority of our samples with reference samples of jemenitica and syriaca haplotypes, 3- jemenitica (Group 2), which includes the far distinct samples (Yem4, Eht5, 4) clustering with haplotypes described to be jemenitica but all sampled from Ethiopia and finally 4-North African haplotypes, mainly characterized by intermissa, Figure 4A. The wide range of genetic diversity recorded for jemenitica populations is indeed unique, even at intra-sample level, Figure 4B. This might not be all surprising for the African lineage which was described to hold the highest genetic diversity among all lineages (Alburaki et al., 2011; Henriques et al., 2022; Franck et al., 2001). Nevertheless, a COI-COII haplotype diversity (H) of 0.980 in a population of 15 samples is unprecedented. The closest (H) that can be reported from previous investigations of African subspecies are 0.908 (Algeria, n = 42) as well as 0.701 (Algeria, n = 582) within intermissa populations (Achou et al., 2015; Ayad et al., 2024), 0.732 (Syria and Lebanon, n = 1,801) within syriaca populations (Alburaki et al., 2011). In contrast, mediocre haplotype diversity (H = 0.597, n = 1,063) was recorded for the honey bee populations of the United States, in which samples overwhelmingly belonged to the North Mediterranean lineage C (Alburaki et al., 2023). While all intermissa and syriaca studies cited above identified some degree of foreign mtDNA C-lineage introgression, our results showed no introgression of any kind to be recorded within the jemenitica populations of the Arabian Peninsula. The absence of introgression was also reported in previous jemenitica investigations from Socotra Island (Yemen) and Saudi Arabia (Alattal et al., 2019; Alattal et al., 2014b). Moreover, the studied populations of jemenitica accepted Tajima’s assumption of neutral evolution (D = - 1.25, p = 0.2), which suggests that these populations are not subjected to systematic human selection, Figure 5. This finding reflects what is observed in this region as most beekeepers practice traditional beekeeping with overwhelmingly no queen importations (Al-Ghamdi and Nuru, 2013b). The haplotype networks are in complete agreement with the COI-COII phylogenetic analysis. This intergenic region indicates that Ethiopian samples (Eth1,4,5) are genetically distant from those of Yemen, Saudi Arabia, Oman, and Jordan, Figure 5. When analyzed with available references, four clusters could be distinguished, Figure 6. These clusters represent two evolutionary lineages (A and M), with a significant intra-A lineage diversity. Two Ethiopian and one Yemeni sample clustered independently with previous reference samples from Ethiopia originating from a central scutellata haplotype (A4). North African intermissa haplotypes (A4p, A9j, A9l) clustered in closer proximity (SNP <4) to the central scutellata (A4) haplotype, Figure 6. The third clade relevant to this study comprises most of the analyzed samples clustering with jemenitica and syriaca reference haplotypes. This cluster suggests that jemenitica populations diverged from the syriaca haplotype (Z8), confirming our previous cyt b finding of intertwining between these two African subspecies despite being geographically distant. To understand this dual relationship, further investigations are needed with a sizable sampling from the contact zone between jemenitica (North) and syriaca (South).

In conclusion, this study provides new insights into the genetic diversity of the Arabian Peninsula’s honey bee populations and two neighboring countries. Despite the limited sample size, our data reveal significant genetic diversity within jemenitica populations. We identified and named eight novel COI-COII haplotypes from Yemen, Saudi Arabia, and Ethiopia that have never been reported. Evidence of subspecies overlapping was recorded, particularly between jemenitica and syriaca. Like many other autochthon subspecies, jemenitica is a remarkable honey bee subspecies that extends on a large geographical area and exemplifies adaptability and resilience in harsh environments. Its ecological importance and roles in both traditional and modern beekeeping underscore the need for continued research and conservation efforts to ensure its survival and support ecosystem health. Finally, our study was based on a limited sample size and focused solely on mtDNA markers. Future research should incorporate genomic DNA analysis with a significantly larger sample size to gain a more comprehensive understanding of the genetic diversity of jemenitica populations.

References

• 1

  AbrahamovichA. H.AtelaO.De La RúaP.GaliánJ. (2007). Assessment of the mitochondrial origin of honey bees from Argentina. J. Apic. Res.46, 191–194. 10.3896/ibra.1.46.3.10

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AchouM.Loucif-AyadL.LegoutH.HmidanH.AlburakiM.GarneryL. (2015). An insightful molecular analysis reveals foreign honeybees among Algerian honeybee populations (Apis mellifera L.). J. Data Min. Genomics and Proteomics6, 1. 10.4172/2153-0602.1000166

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AlattalY.AighamdiA.AlsharhiM. (2014a). Population structure of the Yemeni honey bee (apis mellifera jemenitica) entails an urgent conservation strategy in Saudi Arabia. J. Entomology11, 163–169. 10.3923/je.2014.163.169

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AlattalY.AlsharhiM.AlghamdiA.AlfaifyS.MigdadiH.AnsariM. (2014b). Characterization of the native honey bee subspecies in Saudi Arabia using the mtDNA COI-COII intergenic region and morphometric characteristics. Bull. Insectology67(1), 31–37.

  • Google Scholar
  • View reference in article

• 5

  AlattalY.AlsharhiM.AlghamdiA.FuchsS. (2019). Characterization of Socotra (Yemen) honey bees (Apis mellifera) using morphometric and genetic markers. Bull. Insectology72, 281–285.

  • Google Scholar
  • View reference in article

• 6

  AlburakiM. (2011). Analyse de la diversité moléculaire et morphométrique des populations d'abeilles de Syrie Apis mellifera syriaca: application à la conservation et à la sélection des populations d'abeilles. PhD theses. French Ministry of Higher Education and Research, Sorbonne University. Ph.D.

  • Google Scholar
  • View reference in article

• 7

  AlburakiM.AlburakiA. (2008). Morphometrical study on syrian honeybee (Apis mellifera syriaca). Emir. J. Food Agric.20, 89–93. 10.9755/ejfa.v20i1.5184

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  AlburakiM.BertrandB.LegoutH.MoulinS.AlburakiA.SheppardW. S.et al (2013). A fifth major genetic group among honeybees revealed in Syria. BMC Genet.14, 117. 10.1186/1471-2156-14-117

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  AlburakiM.MadellaS.LopezJ.BougaM.ChenY. P.VanengelsdorpD. (2023). Honey bee populations of the USA display restrictions in their mtDNA haplotype diversity. Front. Genet.13, 1092121. 10.3389/fgene.2022.1092121

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  AlburakiM.MoulinS.LegoutH.AlburakiA.GarneryL. (2011). Mitochondrial structure of eastern honeybee populations from Syria, Lebanon and Iraq. Apidologie42, 628–641. 10.1007/s13592-011-0062-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  AlghamdiA.AlattalY. (2020). The complete mitochondrial genome and phylogenetic analysis of the Arabian Honeybee Apis mellifera jemenitica (Insecta: hymenoptera: Apidae) from Saudi Arabia. Mitochondrial DNA B Resour.5, 3673–3675. 10.1080/23802359.2020.1832598

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  AlghamdiA.AlsharhiM.AlattalY.AdgabaN. (2012). Morphometric diversity of indigenous honeybees, (linnaeus, 1758), in Saudi Arabia (insecta: apidae). Zoology Middle East57, 97–103. 10.1080/09397140.2012.10648968

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  Al-GhamdiA.NuruA. (2013a). Beekeeping in the kingdom of Saudi Arabia opportunities and challenges. Bee World90, 54–57. 10.1080/0005772x.2013.11417543

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  Al-GhamdiA.NuruA. (2013b). Beekeeping in the kingdom of Saudi Arabia past and present practices. Bee World90, 26–29. 10.1080/0005772x.2013.11417527

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  AlqarniA. S.HannanM. A.OwayssA. A.EngelM. S. (2011). The indigenous honey bees of Saudi Arabia (Hymenoptera, Apidae, Apis mellifera jemenitica Ruttner): their natural history and role in beekeeping. Zookeys134, 83–98. 10.3897/zookeys.134.1677

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  AriasM. C.RindererT. E.SheppardW. S. (2006). Further characterization of honey bees from the Iberian Peninsula by allozyme, morphometric and mtDNA haplotype analyses. J. Apic. Res.45, 188–196. 10.3896/ibra.1.45.4.04

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  AriasM. C.SheppardW. S. (1996). Molecular phylogenetics of honey bee subspecies (Apis mellifera L.) inferred from mitochondrial DNA sequence. Mol. Phylogenet Evol.5, 557–566. 10.1006/mpev.1996.0050

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  AviseJ. C. (2000). Phylogeography: the history and formation of species. Cambridge: Harvard University Press, 447.

  • Google Scholar
  • View reference in article

• 19

  AviseJ. C.ArnoldJ.BallR. M.BerminghamE.LambT.NeigelJ. E.et al (1987). Intraspecific phylogeography: the mitochondrial DNA bridge between population genetics and systematics. Annu. Rev. Ecol. Syst.18, 489–522. 10.1146/annurev.ecolsys.18.1.489

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  AyadA. S.HebertM. P. A.DoironJ. A.Loucif-AyadW.DaasT.SmaggheG.et al (2024). Algerian propolis from distinct geographical locations: chemical profiles, antioxidant capacity, cytotoxicity and inhibition of 5-lipoxygenase product biosynthesis. Chem. Biodivers.21, e202301758. 10.1002/cbdv.202301758

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  BoardmanL.EimanifarA.KimballR. T.BraunE. L.FuchsS.GrunewaldB.et al (2020). The complete mitochondrial genome of Apis mellifera jemenitica (Insecta: hymenoptera: Apidae), the Arabian honey bee. Mitochondrial DNA B Resour.5, 875–876. 10.1080/23802359.2020.1717383

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  BougaM.HarizanisP. C.KiliasG.AlahiotisS. (2005). Genetic divergence and phylogenetic relationships of honey bee Apis mellifera (Hymenoptera: apidae) populations from Greece and Cyprus using PCR-RFLP analysis of three mtDNA segments. EDP SCIENCES. 10.1051/apido:2005021

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  ChenC.LiuZ.PanQ.ChenX.WangH.GuoH.et al (2016). Genomic analyses reveal demographic history and temperate adaptation of the newly discovered honey bee subspecies Apis mellifera sinisxinyuan n. ssp. Mol. Biol. Evol.33, 1337–1348. 10.1093/molbev/msw017

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  CornuetJ. M.GarneryL.SolignacM. (1991). Putative origin and function of the intergenic region between COI and COII of Apis mellifera L. mitochondrial DNA. Genetics128, 393–403. 10.1093/genetics/128.2.393

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  CrozierR. H.CrozierY. C. (1992). The cytochrome b and ATPase genes of honeybee mitochondrial DNA. Mol. Biol. Evol.9, 474–482. 10.1093/oxfordjournals.molbev.a040729

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  CrozierY. C.KoulianosS.CrozierR. H. (1991). An improved test for Africanized honeybee mitochondrial DNA. Experientia47, 968–969. 10.1007/BF01929894

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  De La RuaP.SerranoJ.GalianJ. (1998). Mitochondrial DNA variability in the canary islands honeybees (Apis mellifera L.). Mol. Ecol.7, 1543–1547. 10.1046/j.1365-294x.1998.00468.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  De Oliveira FranciscoF.SilvestreD.AriasM. C. (2001). Mitochondrial DNA characterization of five species of Plebeia (Apidae: meliponini): RFLP and restriction maps. Apidologie32, 323–332. 10.1051/apido:2001132

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  DogantzisK. A.ZayedA. (2019). Recent advances in population and quantitative genomics of honey bees. Curr. Opin. Insect Sci.31, 93–98. 10.1016/j.cois.2018.11.010

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  DrummondA. J.AshtonB.BuxtonS.CheungM.CooperA.DuranC.et al (2011). Geneious v5.4. Available online at: http://www.geneious.com/.

  • Google Scholar
  • View reference in article

• 31

  DuttonR. W.RuttnerF.BerkeleyA.ManleyM. J. D. (1981). Observations on the morphology, relationships and ecology of apis-mellifera of Oman. J. Apic. Res.20, 201–214. 10.1080/00218839.1981.11100498

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  EdgarR. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res.32, 1792–1797. 10.1093/nar/gkh340

  • CrossRef
  • Google Scholar

• 33

  El-NiweiriM. A. A.MoritzR. F. A. (2008). Mitochondrial discrimination of honeybees (Apis mellifera) of Sudan. Apidologie39, 566–573. 10.1051/apido:2008039

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  EstoupA.GarneryL.SolignacM.CornuetJ. M. (1995). Microsatellite variation in honey bee (Apis mellifera L.) populations: hierarchical genetic structure and test of the infinite allele and stepwise mutation models. Genetics140, 679–695. 10.1093/genetics/140.2.679

  • CrossRef
  • Google Scholar

• 35

  FranckP.GarneryL.CelebranoG.SolignacM.CornuetJ. M. (2000). Hybrid origins of honeybees from Italy (Apis mellifera ligustica) and sicily (A. m. sicula). Mol. Ecol.9, 907–921. 10.1046/j.1365-294x.2000.00945.x

  • CrossRef
  • Google Scholar

• 36

  FranckP.GarneryL.LoiseauA.OldroydB. P.HepburnH. R.SolignacM.et al (2001). Genetic diversity of the honeybee in Africa: microsatellite and mitochondrial data. Heredity86, 420–430. 10.1046/j.1365-2540.2001.00842.x

  • CrossRef
  • Google Scholar

• 37

  FtayehA.MeixnerM.FuchsS. (1994). Morphometrical investigation in Syrian honeybees. Apidologie25, 396–401. 10.1051/apido:19940406

  • CrossRef
  • Google Scholar

• 38

  GarneryL.CornuetJ. M.SolignacM. (1992). Evolutionary history of the honey bee Apis mellifera inferred from mitochondrial DNA analysis. Mol. Ecol.1, 145–154. 10.1111/j.1365-294x.1992.tb00170.x

  • CrossRef
  • Google Scholar

• 39

  GarneryL.FranckP.BaudryE.VautrinD.CornuetJ. M.SolignacM. (1998a). Genetic diversity of the west European honey bee (Apis mellifera mellifera and A. m. iberica). I. Mitochondrial DNA. Genet. Sel. Evol.30, S31–S47. 10.1186/1297-9686-30-s1-s31

  • CrossRef
  • Google Scholar

• 40

  GarneryL.FranckP.BaudryE.VautrinD.CornuetJ. M.SolignacM. (1998b). Genetic diversity of the west European honey bee (Apis mellifera mellifera and A. m. iberica). II. Microsatellite loci. Genet. Sel. Evol.30, S49–S74. 10.1051/gse:19980703

  • CrossRef
  • Google Scholar

• 41

  GarneryL.MosshineE. H.OldroydB. P.CornuetJ. M. (1995). Mitochondrial-DNA variation in Moroccan and Spanish honey-bee populations. Mol. Ecol.4, 465–472. 10.1111/j.1365-294x.1995.tb00240.x

  • CrossRef
  • Google Scholar

• 42

  GarneryL.SolignacM.CelebranoG.CornuetJ. M. (1993). A simple test using restricted pcr-amplified mitochondrial-DNA to study the genetic-structure of apis-mellifera L. Experientia49, 1016–1021. 10.1007/bf02125651

  • CrossRef
  • Google Scholar

• 43

  GuindonS.DufayardJ.-F.LefortV.AnisimovaM.HordijkW.GascuelO. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol.59, 307–321. 10.1093/sysbio/syq010

  • CrossRef
  • Google Scholar

• 44

  HallH. G.SmithD. R. (1991). Distinguishing African and European honeybee matrilines using amplified mitochondrial DNA. Proc. Natl. Acad. Sci. U. S. A.88, 4548–4552. 10.1073/pnas.88.10.4548

  • CrossRef
  • Google Scholar

• 45

  HarrisonR. G. (1989). Animal mitochondrial DNA as a genetic marker in population and evolutionary biology. Trends Ecol. Evol.4, 6–11. 10.1016/0169-5347(89)90006-2

  • CrossRef
  • Google Scholar

• 46

  HenriquesD.CostaC.RufinoJ.PintoM. A. (2022). The mitochondrial genome of Apis mellifera siciliana. Mitochondrial DNA B Resour.7, 828–830. 10.1080/23802359.2022.2073844

  • CrossRef
  • Google Scholar

• 47

  IlyasovR. A.LeeM.-L.TakahashiJ.-I.KwonH. W.NikolenkoA. G. (2020). A revision of subspecies structure of western honey bee Apis mellifera. Saudi J. Biol. Sci.27, 3615–3621. 10.1016/j.sjbs.2020.08.001

  • CrossRef
  • Google Scholar

• 48

  KandemirI.KenceM.KenceA. (2000). Genetic and morphometric variation in honeybee (Apis mellifera L.) populations of Turkey. Apidologie31, 343–356. 10.1051/apido:2000126

  • CrossRef
  • Google Scholar

• 49

  KandemirI.KenceM.SheppardW. S.KenceA. (2006). Mitochondrial DNA variation in honey bee (Apis mellifera L.) populations from Turkey. J. Apic. Res.45, 33–38. 10.1080/00218839.2006.11101310

  • CrossRef
  • Google Scholar

• 50

  Kauhausen-KellerD.RuttnerF.KellerR. (1997). Morphometric studies on the microtaxonomy of the species Apis mellifera L. Apidologie28, 295–307. 10.1051/apido:19970506

  • CrossRef
  • Google Scholar

• 51

  Loucif-AyadW.AchouM.LegoutH.AlburakiM.GarneryL. (2014). Genetic assessment of Algerian honeybee populations by microsatellite markers. Apidologie46, 392–402. 10.1007/s13592-014-0331-0

  • CrossRef
  • Google Scholar

• 52

  MadellaS.GrubbsK.AlburakiM. (2020). Non-invasive genotyping of honey bee queens Apis mellifera L.: transition of the DraI mtDNA COI-COII test to in silico. Insects12, 19. 10.3390/insects12010019

  • CrossRef
  • Google Scholar

• 53

  MeixnerM. D.AriasM. C.SheppardW. S. (2000). Mitochondrial DNA polymorphisms in honey bee subspecies from Kenya. Apidologie31, 181–190. 10.1051/apido:2000115

  • CrossRef
  • Google Scholar

• 54

  MeixnerM. D.WorobikM.WildeJ.FuchsS.KoenigerN. (2007). Apis mellifera mellifera in eastern Europe †morphometric variation and determination of its range limits. Apidologie38, 191–197. 10.1051/apido:2006068

  • CrossRef
  • Google Scholar

• 55

  MoritzR. F. A.CornuetJ. M.KrygerP.GarneryL.HepburnH. R. (1994). Mitochondrial DNA variability in South African honeybees (Apis mellifera L). Apidologie25, 169–178. 10.1051/apido:19940205

  • CrossRef
  • Google Scholar

• 56

  OldroydB. P.WongsiriS. (2006). Asian honey bees. Biology, conservation and human interactions. Cambridge, Ma: Harvard University Press.

  • Google Scholar

• 57

  PalmerM. R.SmithD. R.KaftanogluO. (2000). Turkish honeybees: genetic variation and evidence for a fourth lineage of Apis mellifera mtDNA. J. Hered.91, 42–46. 10.1093/jhered/91.1.42

  • CrossRef
  • Google Scholar

• 58

  ParadisE. (2010). pegas: an R package for population genetics with an integrated-modular approach. Bioinformatics26, 419–420. 10.1093/bioinformatics/btp696

  • CrossRef
  • Google Scholar

• 59

  ParadisE.SchliepK. (2019). Ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics35, 526–528. 10.1093/bioinformatics/bty633

  • CrossRef
  • Google Scholar

• 60

  PintoM. A.RubinkW. L.CoulsonR. N.PattonJ. C.JohnstonJ. S. (2004). Temporal pattern of africanization in a feral honeybee population from Texas inferred from mitochondrial DNA. Evolution58, 1047–1055. 10.1111/j.0014-3820.2004.tb00438.x

  • CrossRef
  • Google Scholar

• 61

  R CORE TEAM (2022). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/.

  • Google Scholar

• 62

  RindererT. E. (1986). Bee genetics and breeding. Orlando: Academic Press.

  • Google Scholar

• 63

  RortaisA.ArnoldG.AlburakiM.LegoutH.GarneryL. (2011). Review of the DraI COI-COII test for the conservation of the black honeybee (Apis mellifera mellifera). Conserv. Genet. Resour.3, 383–391. 10.1007/s12686-010-9351-x

  • CrossRef
  • Google Scholar

• 64

  RuttnerF. (1988). Biogeography and taxonomy of honeybees. Berlin: Springer-Verlag.

  • Google Scholar

• 65

  RuttnerF.TassencourtL.LouveauxJ. (1978). BIOMETRICAL-STATISTICAL analysis of the geographic variability of APIS mellifera L. I. Material and methods. Apidologie9, 363–381. 10.1051/apido:19780408

  • CrossRef
  • Google Scholar

• 66

  SheppardW. S.MeixnerM. (2003). Apis mellifera pomonella, a new honey bee subspecies from Central Asia. Apidologie34, 367–375. 10.1051/apido:2003037

  • CrossRef
  • Google Scholar

• 67

  SheppardW.AriasM.ShimanukiH. (1994). Determination of mitochondrial DNA haplotypes from sting remnants of the honeybee Apis mellifera (Hymenoptera: apidae). Bull. entomological Res.84, 551–554. 10.1017/s0007485300032806

  • CrossRef
  • Google Scholar

• 68

  SmithD. R.BrownW. M. (1988). Polymorphisms in mitochondrial DNA of European and Africanized honeybees (Apis mellifera). Experientia44, 257–260. 10.1007/BF01941730

  • CrossRef
  • Google Scholar

• 69

  SmithD. R.BrownW. M. (1990). Restriction endonuclease cleavage site and length polymorphisms in mitochondrial DNA of Apis mellifera mellifera and A. M. carnica (hymenoptera: apidae). Ann. Entomol. Soc. Am.83, 81–88. 10.1093/aesa/83.1.81

  • CrossRef
  • Google Scholar

• 70

  SmithD. R.PalopoliM. F.TaylorB. R.GarneryL.CornuetJ. M.SolignacM.et al (1991). Geographical overlap of two mitochondrial genomes in Spanish honeybees (Apis mellifera iberica). J. Hered.82, 96–100. 10.1093/oxfordjournals.jhered.a111062

  • CrossRef
  • Google Scholar

• 71

  SmithD. R.SlaymakerA.PalmerM.KaftanogluO. (1997). Turkish honey bees belong to the east Mediterranean mitochondrial lineage. Apidologie28, 269–274. 10.1051/apido:19970503

  • CrossRef
  • Google Scholar

• 72

  TamuraK.NeiM. (1993). Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol. Evol.10, 512–526. 10.1093/oxfordjournals.molbev.a040023

  • CrossRef
  • Google Scholar

• 73

  TobeS. S.KitchenerA. C.LinacreA. M. T. (2010). Reconstructing mammalian phylogenies: a detailed comparison of the cytochrome b and cytochrome oxidase subunit I mitochondrial genes. PLOS ONE5, e14156. 10.1371/journal.pone.0014156

  • CrossRef
  • Google Scholar

• 74

  ToparslanE.KarabagK.BilgeU. (2020). A workflow with R: phylogenetic analyses and visualizations using mitochondrial cytochrome b gene sequences. PLoS One15, e0243927. 10.1371/journal.pone.0243927

  • CrossRef
  • Google Scholar

• 75

  ZouY.ZhangZ.ZengY.HuH.HaoY.HuangS.et al (2024). Common methods for phylogenetic tree construction and their implementation in R. Bioeng. (Basel)11, 480. 10.3390/bioengineering11050480

  • CrossRef
  • Google Scholar